Pierre Marion
email: pierre.marion@sorbonne-universite.fr

Adeline Fermanian
email: adeline.fermanian@mines-paristech.fr

Gérard Biau
email: gerard.biau@sorbonne-universite.fr

Jean-Philippe Vert
email: jpvert@google.com

Scaling ResNets in the Large-depth Regime

Keywords: ResNets, deep learning theory, neural ODE, neural network initialization, continuous-time models

published or not. The documents may come

Introduction

Deep residual neural networks

Residual neural networks (ResNets), introduced by He et al. (2016a) in the field of computer vision, were the first deep neural network models successfully trained with several thousand layers. Since then, extensive empirical evidence has shown that increasing the depth leads to significant improvements in performance, while raising new challenges in terms of training (e.g., [START_REF] Wang | Deepnet: Scaling transformers to 1,000 layers[END_REF]. From a high-level perspective, the key feature of ResNets is the presence of skip connections between successive layers. In mathematical terms, this means that the (k + 1)-th hidden state h k+1 ∈ R d follows sequentially from the previous hidden state via the recurrence relation

h k+1 = h k + f (h k , θ k+1 ), 0 k L -1, (1)
where f (•, θ k+1 ) is the layer function parameterized by θ k+1 ∈ R p and L is the number of layers. The skip connection corresponds to the addition of h k on the right-hand side of (1), which is absent in classical feedforward networks. This refinement prevents instability issues during training when L is large, provided training is performed in a careful way [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on ImageNet classification[END_REF]. The idea of adding skip connections has become common practice in the field of deep learning, and is today incorporated in many other models such as Transformers in natural language processing (Vaswani et al., 2017). For simplicity, in the rest of the paper, we continue to use the terminology ResNets to denote any architecture of the form (1), keeping in mind that this framework goes beyond the original model of He et al. (2016a).

The most common architectures have 50-150 layers, but ResNets can be trained with depths up to the order of thousand layers (He et al., 2016b). Yet, the training procedure needs to be carefully crafted to avoid vanishing or exploding gradients, particularly as the depth increases. As pointed out by, e.g., Shao et al. (2020), these instabilities are related to a shift in the magnitude of the variance of a signal as it passes through the network. In the original approach of He et al. (2016a), the issue was mitigated by adding a normalization step, called batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF], which rescales the output of each layer via centering and unit variance normalization. However, this normalization stage introduces practical and theoretical difficulties, among which computational overhead and strong dependence on the batch size (see Brock et al., 2021, and the references therein). A widespread alternative to stabilize training in deep models, explored for example by Yang and Schoenholz (2017), Arpit et al. (2019), Zhang et al. (2019b), and De and Smith (2020), is to incorporate a scaling factor α L in front of the residual term in (1), yielding the model

h k+1 = h k + α L f (h k , θ k+1 ), 0 k L -1.
(2)

There is strong evidence that this scaling factor α L should depend on L, without however any consensus to date on the exact form of this dependence, nor on the mathematical grounding of the approach. Thus, despite progresses on the empirical side, the mathematical forces in action behind the stability of deep ResNets are still poorly understood, although they are key to unlock training at arbitrary depth.

Our goal in the present paper is to take a step forward towards a better theoretical understanding of deep ResNets by providing a thorough probabilistic analysis of the sequence (h k ) 0 k L at initialization when L is large, and by leveraging a continuous-time interpretation of model ( 2) via the so-called neural ordinary differential equation (neural ODE, Chen et al., 2018) paradigm. In a nutshell, our results highlight the intimate connection that exists at initialization between stability of the learning process, the regularity of the weights, and the scaling factor α L . We offer in particular a proper mathematical grounding on why and how to choose the parameter α L as a function of the depth L and the distribution of the weights.

Our contributions

Scaling at initialization. The optimal parameters of ResNets are learnt by minimizing some empirical risk function via a gradient descent algorithm. As highlighted for example by Yang and Schoenholz (2017), Hanin and Rolnick (2018), and Arpit et al. (2019), a good parameter initialization of this learning phase plays a major role in the quality of the learnt model, in particular to avoid vanishing gradients and deadlock at initialization, or exploding gradients and quick divergence of the model parameters at the beginning of training. Moreover, a good initialization allows the use of larger learning rates, which have been shown to correlate with better generalization [START_REF] Jastrzkebski | Three factors influencing minima in SGD[END_REF]. It is thus of great interest to study and understand the role played by scaling of deep ResNets at initialization. This is the context in which we place ourselves in the sequel.

At initialization stage, the weights (θ k ) 1 k L are usually chosen as (realizations of) independent and identically distributed (i.i.d.) random variables, which typically follow a uniform or Gaussian distribution on R p . Accordingly, the sequence (h k ) 0 k L that results from the recursion (2) for a given input to the network takes the form of a sequence of random variables that are not i.i.d. but are actually a martingale. Thus, denoting informally by L the differentiable loss associated with the learning task (classification or regression), the distributions of (h k ) 0 k L and ( ∂L ∂h k ) 0 k L as L becomes large carry useful information on the stability of training. For instance, exploding gradients in the backpropagation phase of learning correspond to the fact that, with high probability, ∂L ∂h 0 ∂L ∂h L , where • denotes the Euclidean norm. Our first contribution, in Section 2, is to provide thorough mathematical statements on the behavior of these distributions (both for finite and infinite L), depending on the value of α L . Among other results, we show that only the choice α L ≈ 1 / √ L yields a non-trivial behavior at initialization, thereby confirming empirical findings in the literature (Arpit et al., 2019;De and Smith, 2020). For α L 1 / √ L, the norms explode exponentially fast with L, which is inappropriate for training. For α L 1 / √ L, the network is almost equivalent to identity, that is, h L ≈ h 0 . The analysis of the different cases as a function of α L is mathematically involved and makes extensive use of concentration tools from random matrix theory.

The continuous approach. As noticed by several authors (Chen et al., 2018;[START_REF] Thorpe | Deep limits of residual neural networks[END_REF]E et al., 2019), model (2) with a scaling factor α L = 1 /L (and not 1 / √ L) is formally similar to the discretization of a differential equation. Thus, when L tends to infinity, the weights and hidden states change continuously with the layer according to the equation

dH t dt = f (H t , Θ t ), t ∈ [0, 1]. (3)
Here, time t is the continuous analogue of the layer index k, H : [0, 1] → R d is a continuoustime hidden state, and Θ : [0, 1] → R p a continuous-time parameter. This important connection between ResNets and differential equations has been identified in the past years under the umbrella name of neural ODE. Since the original article of Chen et al. (2018), this point of view has led to the development of a variety of new continuous-time models, together with innovative architectures and efficient training algorithms [START_REF] Chang | AntisymmetricRNN: A dynamical system view on recurrent neural networks[END_REF][START_REF] Grathwohl | FFJORD: Free-form continuous dynamics for scalable reversible generative models[END_REF]Kidger et al., 2021). The neural ODE paradigm also enabled to leverage the rich theory of differential equations to better understand the mechanisms at work behind deep ResNets (E et al., 2019;Fermanian et al., 2021). However, there is a debated question in the neural ODE community about the choice α L = 1 /L, which guarantees convergence of the discrete model (2) to its continuous-time counterpart (3). As a matter of fact, it seems that this choice is guided by more mathematical than practical considerations, and several authors have suggested that it is inconsistent with what is done in practice [START_REF] Cohen | Scaling properties of deep residual networks[END_REF][START_REF] Bayer | Stability of deep neural networks via discrete rough paths[END_REF]. Moreover, letting α L = 1 /L is somewhat contradictory with the results discussed above, which highlighted that the only non-trivial limit at initialization is

α L = 1 / √ L.
Thus, as a second contribution, we clarify the problem in Section 3 by leveraging our previous results on stability. We show that the value α L = 1 / √ L corresponds in the continuous world to a neural stochastic differential equation (SDE) of the form (3), where now Θ : [0, 1] → R p takes the form of a continuous-time stochastic process, typically a Brownian motion. By contrast, we also prove that the neural ODE regime with α L = 1 /L corresponds to the limit of a ResNet, not with i.i.d. weights as considered before, but with more complex and correlated weight distributions. For these weight distributions, the scaling α L = 1 /L is also a critical value between explosion and identity.

Going further, our third contribution is to exhibit in Section 4 a continuous range of regimes that are controlled by the choice of α L (beyond the cases 1 / √ L and 1 /L) and the distribution of (θ k ) 1 k L at initialization, derived from a continuous-time process Θ with a regularity different from a Brownian motion. More precisely, we show experimentally that there is a strong interplay (with the same three cases-explosion, identity mapping, non-trivial behavior) between the choice of α L and the regularity of (θ k ) 1 k L as a function of the layer index k. In addition, empirical evidence suggests that this interplay impacts both the behavior and performance of the networks during training, beyond initialization.

Related work

The choice of scaling for ResNets has been discussed in many papers, without however reaching a clear consensus on the form this scaling factor should take. For instance, Hanin and Rolnick (2018) state that stability requires α L 1 /L, while Zhang et al. (2019b) show that α L 1 / √ L is enough to ensure stability. On the other hand, [START_REF] Cohen | Scaling properties of deep residual networks[END_REF] claim that the scaling factor observed in practice in trained ResNets is of the form 1 /L β with β ≈ 0.7. Other authors have proposed more complex choices for α L (e.g., Zhang et al., 2019a;Shao et al., 2020). Taking another point of view, De and Smith (2020) observe that batch normalization is empirically equivalent to taking a 1 / √ L normalization factor. [START_REF] Bachlechner | ReZero is all you need: Fast convergence at large depth[END_REF] suggest to learn a scaling parameter α k that is allowed to vary from one layer to another, whereas, in (4), α L is kept constant across layers. These authors observe a great acceleration for training compared to traditional ResNets with no scaling. They also suggest a similar architecture for Transformers and then notice that α k ≈ 1 /L at the end of training.

Closest to our analysis at initialization are the papers of Arpit et al. (2019) and Zhang et al. (2019b). Arpit et al. (2019) develop a theoretical analysis based on mean field approximation that suggests that a scaling factor α L = 1 / √ L prevents vanishing/exploding gradients at initialization, and provide experimental evidence that this approach is competitive with batch normalization. However, the authors do not provide rigorous mathematical statements for the three different cases α

L 1 / √ L, α L ≈ 1 / √ L, and α L 1 / √ L,
nor do they highlight the connection with the continuous-time interpretation. Interestingly, the idea of exploiting the martingale structure to analyze the magnitude of the hidden states is present in Zhang et al. (2019b), who study the convergence of gradient descent for over-parameterized ResNets with different values of α L . Nevertheless, they consider a specific model with Gaussian weights, and only provide asymptotic results when both width and depth tend to infinity.

The connection between the choice of scaling and the continuous-time point of view has previously been noticed by Zhang et al. (2019c), then studied in detail by [START_REF] Cohen | Scaling properties of deep residual networks[END_REF]. The latter show that, under assumptions on the form of the weights, it is possible to derive limiting (stochastic or ordinary) differential equations for the hidden states. However, they do not discuss the transition between these two regimes, nor do they link differential equations regimes with the stability of the network.

Scaling at initialization

Our goal in this section is to study the effect of the scaling factor α L on the stability of ResNets at initialization, assuming that the weights are i.i.d. random variables. We start by making more precise the model and the learning problem introduced in (1).

Model and assumptions

Model. The data is a sample of n pairs (x i , y i ) 1 i n , where x i is the input vector in R n in and y i ∈ R nout is the output vector to be predicted. This setting includes regression and classification (after one-hot encoding of the labels). Specifying the informal recurrence (1), for any input x ∈ R n in , we consider the output F π (x) ∈ R nout of the ResNet model defined by

h 0 = Ax, h k+1 = h k + α L V k+1 g(h k , θ k+1 ), 0 k L -1, F π (x) = Bh L , (4)
where α L > 0 is the scaling factor of the ResNet and π = (A, B,

(θ k ) 1 k L , (V k ) 1 k L ) are its parameters, with A ∈ R d×n in , B ∈ R nout×d , θ k ∈ R p and V k ∈ R d×d for k = 1, . . . , L.
The almost-everywhere differentiable function g : R d × R p → R d encodes the choice of architecture. We note that the model includes initial and final linear layers in order to map the input space R n in into the space of hidden states R d , and symmetrically to map the last hidden state h L into the output space R nout . These two transformations are of little interest to us, since we mostly focus on the behavior of the sequence of hidden states (h k ) 0 k L . Let us finally notice that the results of this section can be adapted to hidden layers that do not have the same width, at the cost of increased technicality.

An important feature of model ( 4) is that the layer function takes the form of a matrixvector multiplication, which will prove crucial to make use of concentration results on random matrices. We stress that this setting is standard in practice and that it encompasses many different types of ResNets. It includes for example simple ResNets where g(h, θ) = σ(h) with σ the activation function, and the original ResNets from He et al. (2016a), which have

g(h, θ) = ReLU(W h + b),
where the parameter is a pair θ = (W, b) with W ∈ R d×d a weight matrix and b ∈ R d a bias, and ReLU: x → max(x, 0) is applied element-wise. This setting also includes attention layers, where g corresponds to the scaled dot-product between keys and queries, as well as convolutional layers. Although the assumptions we make later have to be slightly modified to cover this context, the rationale should extend. We leave this extension for future work.

Throughout the article, we let : R nout × R nout → R + be a loss function, differentiable w.r.t. its first parameter, for example the squared loss or the cross-entropy loss. The objective of learning is to find the optimal parameter π that minimizes the empirical risk L (π) = n i=1 (F π (x i ), y i ). Probabilistic setting at initialization. The minimization of the empirical risk is usually performed by stochastic gradient descent or one of its variants (Goodfellow et al., 2016, Chapter 8). The gradient descent is initialized by choosing the weights as (realizations of) i.i.d. random variables. The parameters θ 1 , V 1 , . . . , θ L , V L in model ( 4) are therefore assumed to be an i.i.d. collection of random variables, where we recall that θ k ∈ R p and V k ∈ R d×d parameterize the k-th layer of the network. In this stochastic context, the successive hidden states h 0 , . . . , h L given a fixed input x are also random variables, but their distribution is not i.i.d.-in fact, under our assumptions, this sequence is a martingale. To avoid unnecessary technicalities, we assume that the sequence (h k ) 0 k L is non-atomic. This is for example the case if the distribution of the parameters is absolutely continuous w.r.t. the Lebesgue measure. In particular, this ensures that the sequence (h k ) 0 k L almost surely does not hit the non-differentiability points of g.

It is stressed that the distribution of the parameters are assumed to be independent of the depth, so that all the dependence on L is captured in the scaling factor α L . This model enables us to consider multiple architectures at once, via the function g. By contrast, some authors formulate the problem of scaling as a choice of the variance at initialization (e.g., Yang and Schoenholz, 2017;[START_REF] Wang | Deepnet: Scaling transformers to 1,000 layers[END_REF], which makes the analysis architecturedependent. However, for a given architecture, these two approaches are essentially equivalent since Var(α L V k ) = α 2 L Var(V k ). The quantity h Lh 0 / h 0 carries key information on the behavior of the network at initialization. On the one hand, if h Lh 0 h 0 , the network is essentially equal to the identity function. On the other hand, if h Lh 0 h 0 , the output of the network explodes. An intermediate situation is when h L -h 0 ≈ h 0 . In addition, another source of information is provided by the gradients of the hidden states with respect to the empirical risk L . If ∂L ∂h 0 -∂L ∂h L ∂L ∂h L , the gradients do not change as they flow through the network, which means that the exact same information is backpropagated throughout the network. Conversely, if ∂L ∂h 0 -∂L ∂h L ∂L ∂h L , the gradients explode during backpropagation. By exploiting the martingale structure of ( h k ) 0 k L , as well as state-of-the-art concentration inequalities for random matrices with sub-Gaussian entries, we provide in this section probabilistic bounds on the magnitude of these various quantities.

Assumptions. Some assumptions are needed on the choices of architecture and initialization. Recall that a real-valued random variable X is said to be s 2 sub-Gaussian (van Handel, 2016, Chapter 3) if for all λ ∈ R, E(exp(λX)) exp( λ 2 s 2 /2). The sub-Gaussian property is a constraint on the tail of the probability distribution. As an example, Gaussian random variables on the real line are sub-Gaussian and so are bounded random variables.

The following assumptions will be needed throughout the section: for any 1 k L, (A 1 ) For some s 1, the entries of √ dV k are symmetric i.i.d. s 2 sub-Gaussian random variables, independent of d and L, with unit variance.

(A 2 ) For some C > 0, independent of d and L, and for any h ∈ R d ,

h 2 2 E g(h, θ k ) 2 h 2 and E g(h, θ k ) 8 C h 8 .
Assumption (A 1 ) is mild and satisfied by all initializations used in practice. For example, the classical Glorot initialization [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]-which is the default implementation in the Keras package [START_REF] Chollet | [END_REF]-takes the entries of V k+1 as uniform U(-3/d, 3/d) variables. This means that √ dV k+1 is initialized with U(-√ 3, √ 3) random variables, which satisfy (A 1 ). Other examples include the Gaussian N (0, 1/d) initialization of [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on ImageNet classification[END_REF] and, for example, initialization with Rademacher variables.

The first part of Assumption (A 2 ) ensures that g(•, θ k+1 ) is not too far away from being an isometry in expectation. The second part is more technical and, roughly, allows to upper bound the deviations of the norm of g(h k , θ k+1 ). Our next Proposition 1 shows that most classical ResNet architectures verify Assumption (A 2 ). For the sake of readability, these models, together with their parameters, are summarized in Table 1 below.

Name

Recurrence relation Parameters

res-1 Simple ResNet h k+1 = h k + α L V k+1 σ(h k ) θ k+1 = ∅ res-2 Parametric ResNet h k+1 = h k + α L V k+1 σ(W k+1 h k ) θ k+1 = W k+1 res-3 Classical ResNet h k+1 = h k + α L V k+1 ReLU(W k+1 h k ) θ k+1 = W k+1
Table 1: Examples of ResNet architectures considered in the paper. In the first two cases, the activation function σ is such that, for all x ∈ R,

a|x| |σ(x)| b|x|, 1 / √ 2 a < b 1. In the last two cases, W k+1 ∈ R d×d .
Proposition 1 Let res-1, res-2, and res-3 be the models defined in Table 1. Then (i) Assumption (A 2 ) is satisfied for res-1.

(ii) Assumption (A 2 ) is satisfied for res-2 and res-3, as soon as the entries of √ dW k+1 , 0 k L -1, are symmetric i.i.d. sub-Gaussian random variables, independent of d and L, with unit variance.

In the models res-1 and res-2, σ can be, for instance, taken as the parametric ReLU function, i.e., σ(x) = x + + sx -, where x + (resp. x -) denotes the positive (resp. negative) part and the slope s ∈ [ 1 / √ 2, 1] is a parameter of the model. Observe also that res-2 differs from res-3 since the classical ReLU function is defined by ReLU(x) = x + and thus does not satisfy the condition |σ(x)| a|x|. Note that there is no bias term in these three models, as this term is commonly initialized to zero, and we are interested in the behavior at initialization.

Probabilistic bounds on the norm of the hidden states

The next two propositions describe how the quantity h Lh 0 / h 0 changes as a function of Lα 2 L . Proposition 2 provides a high-probability bound of interest when Lα 2 L 1. In this case, we see that, with high probability, the network acts as the identity function, directly mapping h 0 to h L . On the other hand, Proposition 3 provides information in the two cases Lα 2 L 1 and Lα 2 L ≈ 1. When Lα 2 L 1, the lower bound (i) indicates an explosion with high probability of the norm of the last hidden state. On the other hand, when Lα 2 L ≈ 1, the bounds (i) and (ii) show that h L randomly varies around h 0 with fluctuation sizes bounded from below and above. Proposition 2 Consider a ResNet (4) such that Assumptions (A 1 ) and(A 2

) are satisfied. If Lα 2 L 1, then, for any δ ∈ (0, 1), with probability at least 1 -δ, h L -h 0 2 h 0 2 2Lα 2 L δ .
Proposition 3 Consider a ResNet (4) such that Assumptions (A 1 ) and (A 2 ) are satisfied.

(i) Assume that d 64 and α 2 L 2 ( √ Cs 4 +4 √ C+16s 4 )d . Then, for any δ ∈ (0, 1), with probability at least 1 -δ, h L -h 0 2 h 0 2 > exp 3Lα 2 L 8 - 11Lα 2 L dδ -1, provided that 2L exp - d 64α 2 L s 2 δ 11 . (5)
(ii) Assume that α 2 L 1 √ C(d+128s 4 ) . Then, for any δ ∈ (0, 1), with probability at least 1δ,

h L -h 0 2 h 0 2 < exp Lα 2 L + 5Lα 2 L dδ + 1.
Note that the assumptions of Proposition 3 on d and α L are mild, since in the learning tasks where deep ResNets are involved, one typically has α L = 1 /L β with β > 0, d 10 2 and L 10 2 . Note also that condition (5) is not severe since, when d and L are large, it encompasses all reasonable values of δ. Propositions 2 and 3 are interesting in the sense that they provide finite-depth high-probability bounds on the behavior of the hidden states, depending on the magnitude of Lα 2 L . The results become clearer by letting α L = 1 /L β , with β > 0, as shown in the following corollary.

Corollary 4 Consider a ResNet (4) such that Assumptions (A 1 ) and (A 2 ) are satisfied, and let (iii

α L = 1 /L β , with β > 0. (i) If β > 1 /2, then h L -h 0 h 0 P ----→ L→∞ 0. (ii) If β < 1 /2 and d 9, then h L -h 0 h 0 P ----→ L→∞ ∞.
) If β = 1 /2, d 64, L ( 1 2 √ Cs 4 + 2 √ C + 8s 4 )d + 96 √ Cs 4
, then, for any δ ∈ (0, 1), with probability at least 1δ,

exp 3 8 - 22 dδ -1 < h L -h 0 2 h 0 2 < exp 1 + 10 dδ + 1, provided that 2L exp - Ld 64s 2 δ 11 .
Corollary 4 highlights three different asymptotic behaviors for h L , depending on the values of β. For β > 1 /2, statement (i) tells that h L converges towards h 0 in probability, as L tends to infinity, which means that the neural network is essentially equivalent to an identity mapping. On the other hand, for β < 1 /2, the norm of h L explodes with high probability. Finally, for the critical value β = 1 /2, we see that h L fluctuates around h 0 , with a fluctuation size independent of L. Observe that the lower bound in (iii) is not trivial as soon as exp( 3 /8 -11 /dδ) > 1, i.e., d > 99 /64δ. The message of Corollary 4 is that the only scaling leading to a non-degenerate distribution at initialization is for β = 1 /2.

The three statements of Corollary 4 are illustrated in Figure 1. In this experiment, we consider model res-3, a random Gaussian observation x in dimension n in = 64, and parameters initialized with a uniform distribution U(-3/d, 3/d). We refer to Appendix E for a detailed setup of all the experiments of the paper. Figure 2a shows the empirical distribution of h L / h 0 when β = 1 /2 for a large number of realizations. This figure illustrates in particular that our bounds are reasonably sharp, since the bounds indicate that the first quartile of the distribution is larger than 0.87 (whereas the first quartile of the empirical histogram is equal to 1.21) and the third quartile is less than 2.06 (whereas the third quartile of the empirical histogram is equal to 1.34). Determining the exact distribution of h L / h 0 is an interesting avenue for future research that is beyond the scope of the present article. There is however a strong indication that the ratio follows a log-normal distribution, as confirmed by a normality test on (the log of) the empirical distribution. In a nutshell, the proofs of Propositions 2 and 3 rest upon controlling of the norm of the hidden states, which obeys the recurrence

h k+1 2 = h k 2 + α 2 L V k+1 g(h k , θ k+1 ) 2 + 2α L h k , V k+1 g(h k , θ k+1 ) , (6)
where •, • denotes the standard scalar product in R d . Taking the expectations on both side, one deduces with Assumptions (A 1 ) and (A 2 ) that

E V k+1 g(h k , θ k+1 ) 2 = E E V k+1 g(h k , θ k+1 ) 2 h k , θ k+1 = E g(h k , θ k+1 ) 2 ≈ h k 2 (7)
and

E h k , V k+1 g(h k , θ k+1 ) = E E h k , V k+1 g(h k , θ k+1 ) h k , θ k+1 = 0. ( 8
)
The equalities ( 7) and ( 8) allow deriving without further work bounds in expectation on h L , as already observed by Arpit et al. (2019). However, the results we are after are stronger since they involve high-probability bounds. A finer control of the deviations of V k+1 g(h k , θ k+1 ) 2 and h k , V k+1 g(h k , θ k+1 ) is then needed. This involves concentration inequalities on random matrices with sub-Gaussian entries.

Probabilistic bounds on the gradients

Propositions 2 and 3 provide insights on the output of the network when L is large. However, they do not carry information on the backwards dynamics of propagation of the gradients of the loss p k = ∂L ∂h k ∈ R d . Assessing the dynamics of the (p k ) 0 k L as a function of L is important since the behavior of this sequence impacts trainability of the network at initialization. Thus, in this subsection, we are interested in quantifying the magnitude of p 0p L / p L , when L is large. Notice that, contrarily to the previous subsection where we were mostly interested in the last hidden state h L , the quantity of interest is now p 0 (not p L ), the gradient at index 0. The reason is that the sequence (p k ) 0 k L is defined backwardly, as we will see below. We also stress that (p k ) 0 k L is the sequence of derivatives of the loss w.r.t. the hidden states h k , and not w.r.t. the parameters. The reason for considering this sequence is that the p k are involved in the backpropagation algorithm and are therefore essential for assessing the stability of the gradient descent (see, e.g., Arpit et al., 2019).

Analyzing the behavior of the sequence (p k ) 0 k L is challenging since, according to the backpropagation (or reverse-mode differentiation) formula, one has

p k = p k+1 + α L ∂g(h k , θ k+1 ) ∂h V k+1 p k+1 .
Taking the norm,

p k 2 = p k+1 2 + α 2 L ∂g(h k , θ k+1 ) ∂h V k+1 p k+1 2 + 2α L p k+1 , ∂g(h k , θ k+1 ) ∂h V k+1 p k+1 .
Although the equation looks qualitatively similar to ( 6), it has the unpleasant feature that

∂g(h k ,θ k+1 ) ∂h depends on h k , hence on θ 1 , V 1 , . . . , θ k , V k , while p k+1 depends on θ k+2 , V k+2 , . . . , θ L-1 , V L-1
. This forbids applying directly the same proof techniques as for the hidden states. Therefore, to extract useful information from this recurrence equation, one needs to characterize the dependence of the distribution of

∂g(h k ,θ k+1 ) ∂h
with respect to h k . To do so, it is sometimes assumed that these two quantities are independent (see, e.g., Yang and Schoenholz, 2017). However, assuming independence remains a strong requirement, which is not verified for many network architectures (for example model res-1). We tackle the problem from a different point of view and propose an alternative approach based on forward-mode differentiation, valid under a much weaker assumption. The cost we pay is that we obtain results in expectation and not in high probability.

Let us sketch our approach before stating the results more formally. Let z ∈ R d , and, for any 0 i, j L, let ∂h j ∂h i ∈ R d×d be the Jacobian matrix of h j with respect to h i . Recall that the (m, n)-th entry of this matrix equals the derivative of the m-th coordinate of h j w.r.t. the n-th coordinate of h i . Then, letting q k (z) = ∂h k ∂h 0 z, we have, by the chain rule,

q k+1 (z) = ∂h k+1 ∂h k q k (z) = q k (z) + α L V k+1 ∂g(h k , θ k+1 ) ∂h q k (z). (9)
Identity ( 9), which is similar to (4), expresses q k+1 (z) as a function of q k (z), and therefore respects the flow of information. Next, assuming that z is random with a Gaussian distribution, it is possible to express one of our quantities of interest, p 0 / p L , as a function of the last vector q L (z). Indeed,

p 0 2 p L 2 = 1 p L 2 E z∼N (0,I d ) p 0 z 2 = E z∼N (0,I d ) p L p L q L (z) 2 , ( 10
)
where I d is the identity matrix in R d and the second equality is a consequence of

p 0 z = ∂L ∂h 0 z = ∂L ∂h L ∂h L ∂h 0 z = p L q L (z).
In summary, the recurrence (9) allows us to derive bounds on the norm of q L (z), which can then transfer to p 0 / p L via (10). For this, it is necessary to make the following assumption on the ratio p L / p L :

(A 3 ) Let b = p L / p L . Then E(b|h L ) = 0 and E(b b|h L ) = I d /d.
It is a mild assumption, which is verified for instance if n out = 1 with squared error (for regression) or cross-entropy (for binary classification). In these cases, p L / p L = B / B F , where • F is the Frobenius norm and B is the weight matrix of the last layer. We finally need the following assumption, which is the equivalent of Assumption (A 2 ) for the gradients.

(A 4 ) One has, almost surely,

q k 2 2 E ∂g(h k , θ k+1 ) ∂h q k 2 h k , q k q k 2 .
Assumption (A 4 ) is satisfied by all the standard architectures listed in Table 1, as shown by the next proposition.

Proposition 5 Let res-1, res-2, and res-3 be the models defined in Table 1. Assume that (A 1 ) is satisfied and σ is almost everywhere differentiable, with a σ b. Then (i) Assumption (A 4 ) is satisfied for res-1.

(ii) Assumption (A 4 ) is satisfied for res-2 and res-3, when the entries of √ dW k+1 , 0 k L -1, are symmetric i.i.d. random variables, independent of d and L, with unit variance.

The next two propositions are the counterparts of Proposition 2 and Proposition 3 for the gradient dynamics.

Proposition 6 Consider a ResNet (4) such that Assumptions (A 1 )-(A 4 ) are satisfied. If Lα 2 L 1, then, for any δ ∈ (0, 1), with probability at least 1δ,

p 0 -p L 2 p L 2 2Lα 2 L δ .
Proposition 7 Consider a ResNet (4) such that Assumptions (A 1 )-(A 4 ) are satisfied. Then

1 + 1 2 α 2 L L -1 E p 0 -p L 2 p L 2 (1 + α 2 L ) L -1.
A simple corollary of the propositions above is as follows.

Corollary 8 Consider a ResNet (4) such that Assumptions (A 1 )-(A 4 ) are satisfied, and take α L = 1 /L β , with β > 0. Then (ii

(i) If β > 1 /2, p 0 -p L p L P ----→ L→∞ 0. (ii) If β < 1 /2, E p 0 -p L 2 p L 2 L→∞ ----→ ∞.
) If β = 1 /2, exp 1 2 -1 E p 0 -p L 2 p L 2 exp(4) -1.
Corollary 8 is illustrated in Figure 3. The experimental protocol is the same as in Figure 1, but we now track p 0 and p L , the gradients of the loss L with respect to the first and the last hidden states. In accordance with our results, when β > 1 /2, the gradient remains the same from one layer to another (left plot). On the other hand, the middle plot clearly shows that when β < 1 /2 the gradient explodes. Once again, the case β = 1 /2 (right plot) is the only one for which the distribution of gradients at initialization is non-trivial. Figure 2b illustrates that the empirical distribution of gradients in this case also seems to be log-normal.

In summary, this and the previous subsection both point towards the same conclusion: there are three different cases, depending on the value of β-explosion when β < 1 /2, nondegenerate limit when β = 1 /2, and identity when β > 1 /2. In the explosion case, it is well known that the network cannot be trained (Yang and Schoenholz, 2017). The theory thus points out that the value 1 /2 plays a pivotal role. Remarkably, this value has a specific interpretation in the continuous-time point of view of ResNets, in terms of SDE. This is the topic that we address in the next section.

Scaling in the continuous-time setting

Starting with the discrete ResNet (4), it is tempting to let L go to infinity and consider the network as the discretization of a differential equation where the layer index k ∈ {0, . . . , L} is replaced by the time index t ∈ [0, 1]. This interpretation of deep neural networks has been popularized by Chen et al. (2018) and is often referred to as the neural ODE paradigm. Notice that this setting is different from the so-called mean-field analysis, where the width of the network is assumed to be infinite beforehand. In our setting, the width d is assumed to be finite and fixed.

Convergence towards an SDE in the large-depth regime

One of the main messages of Section 2 is that the standard initialization with i.i.d. parameters leads to a non-degenerate model for large values of L only if Lα 2 L ≈ 1 (Propositions 2 and 3), or, equivalently, if β = 1 /2 when α L = 1 /L β (Corollary 4). Remarkably, in the continuous-time limit, this regime corresponds to the discretization of an SDE. Indeed, consider for simplicity the (discrete) ResNet model res-1

h k+1 = h k + 1 √ L V k+1 σ(h k ), 0 k L -1, ( 11
)
where the entries of V k+1 are assumed to be i.i.d. N (0, 2 /d). Recall the following definition:

Definition 9 A one-dimensional Brownian motion (B t ) t∈[0,1] is a continuous-time stochastic process with B 0 = 0, almost surely continuous, with independent increments, and such that for any 0 s < t 1, B t -B s ∼ N (0, ts).

Now, let B : [0, 1] → R d×d be a (d × d)-dimensional Brownian motion, in the sense that the (B ij ) 1 i,j d are independent one-dimensional Brownian motions. Thus, for any 0 k L -1 and any 1 i, j d, we have

B(k+1) /L,i,j -Bk /L,i,j ∼ N 0, 1 L ,
and the increments for different values of (i, j, k) are independent. As a consequence, the recurrence ( 11) is equivalent in distribution to the recurrence

h k+1 = h k + 2 d σ(h k )(B(k+1) /L -Bk /L ), 0 k L -1.
(Note that this is true because V k+1 has the same distribution as V k+1 .) We recognize the Euler-Maruyama discretization [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] on the {k/L, 0 k L} mesh of the SDE

dH t = 2 d σ(H t )dB t , t ∈ [0, 1], ( 12
)
where the output of the network is now a function of the final value of H, that is, H 1 . The link between the discrete ResNet (11) and the SDE ( 12) is formalized in the next proposition.

Proposition 10 Consider the res-1 model, where the entries of V k+1 are i.i.d. Gaussian N (0, 2 /d) random variables. Assume that the activation function σ is Lipschitz continuous.

Then the SDE (12) has a unique solution H and, for any 0 k L,

E Hk /L -h k c √ L ,
for some c > 0.

Notice that the requirement that σ is Lipschitz continuous is satisfied by most classical activation functions, including ReLU. This proposition is interesting for several reasons. First, the scaling β = 1 /2, which is exactly the one that yields a non-trivial dynamics at initialization, corresponds in the continuous world to a remarkably 'simple' model of diffusion. This shows that very deep neural networks properly initialized with i.i.d. weights are equivalent to solutions of SDE. This analogy opens interesting perspectives for training deep networks using automatic differentiation for solutions of neural SDE [START_REF] Li | Scalable gradients and variational inference for stochastic differential equations[END_REF].

Second, we stress that the emergence of an SDE instead of an ODE carries an important message. Several authors (including, e.g., [START_REF] Thorpe | Deep limits of residual neural networks[END_REF] have shown that, under appropriate assumptions, a deep ResNet converges in the large depth limit to an ODE and not an SDE. The reason why we obtain an SDE here is intrinsically connected with the choice of i.i.d. initialization for the weights, which makes a Brownian motion appear at the limit, as highlighted above. In other words, the i.i.d. initialization, the choice β = 1 /2 (the relevant critical value exhibited in Section 2), and the emergence of an SDE are intimately linked together. On the other hand, the case β = 1 matches with an ODE if the initialization is not i.i.d., as we will see in Subsection 3.2.

Finally, we point out that Proposition 10 states the convergence of a ResNet towards an SDE for the basic architecture res-1 and for Gaussian initialization. The extension to more general settings is an interesting direction of research, although clearly beyond the scope of the present paper (see, e.g., [START_REF] Peluchetti | Infinitely deep neural networks as diffusion processes[END_REF][START_REF] Peluchetti | Infinitely deep neural networks as diffusion processes[END_REF][START_REF] Cohen | Scaling properties of deep residual networks[END_REF], for results in this direction).

Scaling in the neural ODE setting

Convergence towards an ODE. The basic message of our Proposition 10 is that an i.i.d. initialization, together with β = 1 /2, leads to an SDE rather than an ODE. A natural question is then whether a different choice of weight distributions (at initialization) and scaling can lead to a classical neural ODE.

To answer this question and leave the world of i.i.d. initialization, we assume that the weights (V k ) 1 k L and (θ k ) 1 k L are discretizations of smooth functions V : [0, 1] → R d×d and Θ : [0, 1] → R p . We then consider the general iteration (4) with α L = 1 /L, that is,

h 0 = Ax, h k+1 = h k + 1 L V k+1 g(h k , θ k+1 ), 0 k L -1, (13)
where

V k = V k/L and θ k = Θ k/L . Of course, it is still possible to consider (V k ) 1 k L (resp. (θ k ) 1 k L ) as random variables, by letting (V t ) t∈[0,1] (resp. (Θ t ) t∈[0,1]
) be a continuous-time stochastic process. In this model, we shall need the following assumption:

(A 5 ) For any 0 k L -1, one has V k = V k/L and θ k = Θ k/L
, where the stochastic processes V and Θ are almost surely Lipschitz continuous and bounded.

More precisely, almost surely, there exist K

V , K Θ , C V , C Θ > 0, such that, for any s, t ∈ [0, 1], V t -V s K V |t -s|, Θ t -Θ s K Θ |t -s|, V t C V , Θ t C Θ .
A typical model that satisfies Assumption (A 5 ) is obtained by letting the entries of V and Θ be independent Gaussian processes with expectation zero and squared exponential covariance K(x, x ) = exp(-(x-x ) 2 2 2 ), where > 0.

We shall also need the following requirement on g, which is satisfied by all our models as soon as σ is Lipschitz continuous: (A 6 ) The function g is Lipschitz continuous on compact sets, in the sense that for any compact P ⊆ R p , there exists K P > 0 such that, for all h, h ∈ R d , θ ∈ P,

g(h, θ) -g(h , θ) K P h -h ,
and for any compact D ⊆ R d , there exists K D,P > 0 such that, for all h ∈ D, θ, θ ∈ P,

g(h, θ) -g(h, θ ) K D,P θ -θ .
Under Assumptions (A 5 ) and (A 6 ), the recurrence ( 13) almost surely converges towards the neural ODE given by

H 0 = Ax, dH t = V t g(H t , Θ t )dt, t ∈ [0, 1], ( 14
)
as shown by the proposition below.

Proposition 11 Consider model (13) such that Assumptions (A 5 ) and (A 6 ) are satisfied.

Then the ODE (14) has a unique solution H, and, almost surely, there exists some c > 0 such that, for any

0 k L, Hk /L -h k c L .
It should be stressed that the transition from the discrete recurrence (13) to the continuous-time differential equation ( 14) relies on the assumptions that the weight sequences (θ k ) 1 k L and (V k ) 1 k L are the discretizations of smooth limiting processes Θ and V on the one hand, and that the scaling α L is chosen as 1 /L on the other hand. From a practical perspective, Proposition 11 shows that it is possible to initialize ResNets in the ODE regime, by choosing a smooth stochastic process, discretizing it at each layer, and taking a 1 /L scaling. This is in sharp contrast with the results of Sections 2 and 3.1, which show that the usual i.i.d. procedure leads to a neural SDE.

Stability and scaling. Assuming that the weights of the network are discretizations of a smooth function (Assumption (A 5 )), it is possible to obtain stability results, depending on the value of β, similarly to what has been done in Section 2. We show below that β = 1 is a critical value, by examining the hidden states, in the same way as β = 1 /2 is a critical value in the i.i.d. setting. Similar results can be shown for the gradients. We begin by a proposition handling the cases β > 1 and β = 1.

Proposition 12 Consider a ResNet (4) such that Assumptions (A 5 ) and (A 6 ) are satisfied. (ii) If β = 1, then, almost surely, there exists some c > 0 such that

Let α L = 1 /L β , with β > 0. (i) If β > 1, then, almost surely, h L -h 0 h 0 L→∞ ----→ 0.
h L -h 0 h 0 c.
The explosion case (β < 1) is more delicate to deal with. We prove it for a linear model, and leave for future work the extension to more general cases.

Proposition 13 Consider the res-1 model, taking σ as the identity function. Assume that Assumption (A 5 ) is satisfied and that V T 0 has a positive eigenvalue. Let α L = 1 /L β , with β ∈ (0, 1). Then, almost surely,

max k h k -h 0 h 0 L→∞ ----→ ∞.
The assumption of the existence of a positive eigenvalue for V 0 is mild. For instance, if the entries of V 0 are i.i.d. random variables with finite moments of all order, [START_REF] Götze | Rate of convergence to the circular law via smoothing inequalities for log-potentials[END_REF] show that such an eigenvalue exists with probability at least 1 -1 /d for d large enough.

In this setting, we observe experimentally a behavior of the output and of the gradients when L grows large similar to the one explored in Section 2. This is illustrated in Figures 4 and5, which mirror Figures 1 and3 in Section 2. The figures clearly show that there exist three cases for the output and for the gradients: an identity case (left plots), an explosion case (middle), and a non-trivial case separating explosion and identity (right). However, the remarkable point is that the separation occurs for β = 1, and not β = 1 /2, as predicted by Propositions 12 and 13.

Experiments

We experimentally investigate in this section two questions. The first one is to know whether there exists a range of scaling factors β > 0 and weight initializations, beyond the i.i.d. and the smooth regimes. The second question is whether our analysis, which pertains to the initialization phase, provides insights into the training phase, beyond initialization.

Intermediate regimes

In order to describe the transition between the i.i.d. and smooth cases, a possible route is to consider that the weights are increments of a γ-Hölder stochastic process. This model is interesting insofar as the Brownian motion (SDE regime) is ( 1 /2ε)-Hölder (ε > 0) and a Lipschitz continuous stochastic process (ODE regime) is 1-Hölder.

In line with the above, in a series of experiments, we initialize the weights as increments of a fractional Brownian motion (B H t ) t∈[0,1] . Recall that B H is a continuous-time Gaussian process, starting at zero, with zero expectation for all t ∈ [0, 1], and covariance function

E(B H s B H t ) = 1 2 (|s| 2H + |t| 2H -|t -s| 2H ), 0 s, t 1,
where H ∈ (0, 1) is called the Hurst index. This index describes the raggedness of the process, with a higher value leading to a smoother process. When H = 1 /2, the process is a standard Brownian motion (Definition 9), whose increments are independent by construction. When H > 1 /2, the increments of the process are positively correlated, while if H < 1 /2 the increments are negatively correlated. Importantly, a fractional Brownian motion with Hurst index H is (Hε)-Hölder continuous for any ε > 0. In the limit when H → 1, the trajectories converge to linear functions (whose increments satisfy (A 5 )). As an illustration, Figure 6 depicts three realizations of a fractional Brownian motion with H = 0.2 (left), H = 0.5 (middle), and H = 0.8 (right).

In order to assess the effect of the scaling factor β and the Hurst index H, we initialize a neural network res-3 with d = 40, L = 1000, various values of β ∈ [0.2, 1.3], and with weights taken as increments of fractional Brownian motions with various Hurst indices H ∈ (0, 1). Figure 7 depicts the empirical magnitude of the output and the gradients at initialization as a function of the Hurst index H and the scaling factor β. First note that we Figure 7: Magnitude of the outputs and of the gradients as a function of the regularity of the weights (Hurst index H) and of the scaling factor β. The orange zone corresponds to the explosion regime, i.e., h Lh 0 h 0 and p 0p L p L . The blue zone corresponds to the identity regime, i.e., h Lh 0 h 0 and p 0p L p L . Finally, the black zone is an intermediate regime, where h Lh 0 ≈ h 0 and p 0p L ≈ p L . recover the two regimes (i.i.d. and smooth) discussed so far. For H = 1 /2, the i.i.d. regime kicks in, with explosion (β < 1 /2, orange zone), non-trivial behavior (β = 1 /2, black zone), and identity (β > 1 /2, blue zone). Likewise, we see at H = 1 a similar pattern in the smooth regime, with, as predicted by Proposition 12, a critical value β = 1. Beyond these two specific cases, we observe for an index H varying in ( 1 /2, 1) a whole range of intermediate situations, where the transition between identity and explosion seems to happen for a critical β = H. Interestingly, for H < 1 /2, the transition seems to saturate at the value β = 1 /2.

The take-home message is that the choice of the scaling of a ResNet seems to be closely linked to the regularity of the weights as a function of the layer. More precisely, for all regimes, the critical scaling factor between explosion and identity seems to have a natural interpretation as the (Hölder) regularity of the underlying continuous-time stochastic process.

We believe that the mathematical understanding of this connection, beyond the fractional Brownian motion case, is a promising research direction for the future. Finally, these experiments suggest that it is sensible to initialize a ResNet for any value of the scaling β ∈ ( 1 /2, 1), while avoiding the identity and explosion situations, by simulating a fractional Brownian motion of Hurst index H = β and initializing the weights as the increments of this process.

Beyond initialization

At initialization, before the gradient descent, the distribution of the weights (θ k ) 1 k L and (V k ) 1 k L is chosen by the practitioner. By contrast, during and after training, control is lost on these distributions, making the picture more complex. In particular, the existence and characterization of a continuous-time stochastic process whose discretization matches the trained ResNet is an interesting but difficult problem. Attacking this question requires a fine understanding of the interaction between training dynamics and the regularity of the sequence of the weights during the gradient descent. However, there is experimental evidence that the trained weights exhibit strong structure as a function of the layer index k [START_REF] Cohen | Scaling properties of deep residual networks[END_REF][START_REF] Bayer | Stability of deep neural networks via discrete rough paths[END_REF], and that their regularity strongly depends on the choice of initialization. Figure 8 depicts this mechanism by plotting a given coordinate of θ k as a function of the layer index k ranging from 1 to the depth L = 1000, after training. To investigate the link between regularity of the weights at initialization, scaling, and performance after training, we train ResNets on the datasets MNIST [START_REF] Deng | The MNIST database of handwritten digit images for machine learning research[END_REF] and CIFAR-10 (Krizhevsky, 2009). As in Subsection 4.1, we initialize the ResNets with various scaling factors and weights that are increments of fractional Brownian motions with different regularities. Then, for each combination of weight initialization and scaling factor, the ResNet is trained using the Adam optimizer [START_REF] Kingma | A method for stochastic optimization[END_REF] for 10 epochs. The results in terms of accuracy are presented in Figure 9 (light orange = good performance, blue = bad performance). We observe a pattern similar to the one of Figure 7, however shifted downwards. This means that, for a given regularity, the network is unable to learn if it is initialized with a scaling too far below the critical value, which of course is connected with the gradient explosion issue discussed previously. On the other hand, and perhaps more surprisingly, the performance seems to be more or less stable in the identity region, with perhaps a small degradation in the case of CIFAR-10. This somewhat contrasts with the results from Yang and Schoenholz (2017), who exhibit a decrease in performance for i.i.d. initialization and a large scaling factor β. Note however that, in our experiments, we adapted the learning rate of the gradient descent on a grid by cross-validation. This was done to prevent a slowdown in training when the scaling factor β is too large, since, in the gradient descent, the gradients updates are also scaled by the factor 1 /L β -which gets smaller as β increases. The interplay between the learning rate and the scaling factor is one of the keys to better assess how the performance of the trained network is connected with the scaling. With respect to statement (ii), it is enough to show that for any h ∈ R d and any random matrix W satisfying the assumptions of the proposition, one has

h 2 2 E σ(W h) 2 h 2 and E σ(W h) 8 C h 8 ,
as well as

E ReLU(W h) 2 = h 2 2 and E ReLU(W h) 8 C h 8 .
The two claims with the squared norms are consequences of Lemmas 16 and 17 in Appendix B, together with the fact that the variance of the entries of W equals 1/d. In order to prove the other two statements, first note that

E( σ(W h) 8 ) E( W h 8 ) and E( ReLU(W h) 8 ) E( W h 8 ).
The results are then consequences of Lemma 21 in Appendix C, which states that

E W h 8 h 8 1 + 384s 4 d + 3072s 6 d 2 1 + 384s 4 + 3072s 6 .

A.2 Proof of Proposition 2

According to Lemma 14 below, one has

E h L -h 0 2 h 0 2 1 + α 2 L L -1 . But, for Lα 2 L 1, we have (1 + α 2 L ) L -1 exp(Lα 2 L ) -1 2Lα 2 L . Therefore, E h L -h 0 2 h 0 2 2Lα 2 L ,
and the result follows from Markov's inequality.

Lemma 14 Consider a ResNet (4) such that Assumptions (A 1 ) and (A 2 ) are satisfied.

1 + α 2 L 2 L -1 E h L -h 0 2 h 0 2 1 + α 2 L L -1 .
Proof (Lemma 14) Taking the squared norm of the forward update rule (4) and dividing by h 0 2 yields

h k+1 2 h 0 2 = 1 h 0 2 h k 2 + α 2 L V k+1 g(h k , θ k+1 ) 2 + 2α L h k , V k+1 g(h k , θ k+1 ) . (15)
We deduce by Assumptions (A 1 ) and (A 2 ) that

1 + α 2 L 2 E h k 2 h 0 2 E h k+1 2 h 0 2 (1 + α 2 L )E h k 2 h 0 2 .
Therefore, by recurrence, we are led to

1 + α 2 L 2 k E h k 2 h 0 2 (1 + α 2 L ) k . (16) Now, observe that h L = h 0 + α L L-1 k=0 V k+1 g(h k , θ k+1
). Thus, we have

E h L -h 0 2 h 0 2 = α 2 L L-1 k,k =0 E g(h k , θ k+1 ) V k+1 V k +1 g(h k , θ k +1 ) h 0 2 .
By conditioning on all random variables except V k +1 for k < k (and V k+1 for k > k ), it is easy to see that the only non-zero terms are when k = k . This yields

E h L -h 0 2 h 0 2 = α 2 L L-1 k=0 E V k+1 g(h k , θ k+1 ) 2 h 0 2 α 2 L L-1 k=0 E h k 2 h 0 2
(by Assumptions A 1 and A 2 )

α 2 L L-1 k=0 (1 + α 2 L ) k
(by ( 16))

= 1 + α 2 L L -1 .
Similarly,

E h L -h 0 2 h 0 2 α 2 L 2 L-1 k=0 E h k 2 h 0 2 = 1 + α 2 L 2 L -1 .

A.3 Proof of Proposition 3

Dividing (15) by h k 2 and taking the logarithm leads to

ln( h k+1 2 ) = ln( h k 2 ) + ln 1 + α 2 L V k+1 g(h k , θ k+1 ) 2 h k 2 + 2α L h k h k , V k+1 g(h k , θ k+1 ) h k . Let Y k,1 = α 2 L V k+1 g(h k , θ k+1 ) 2 h k 2 , Y k,2 = 2α L h k h k , V k+1 g(h k , θ k+1 ) h k ,
and

Y k = Y k,1 + Y k,2
. The proof of Proposition 3 strongly relies on the following lemma, which provides technical information on the moments of Y k,1 and Y k,2 . For the sake of clarity, its proof is postponed to Appendix B.

Lemma 15 Assume that Assumptions (A 1 ) and (A 2 ) are satisfied. Then

(E 1 ) E(Y k,2 |h k ) = 0. (E 2 ) α 2 L 2 E(Y k,1 |h k ) α 2 L . (E 3 ) E(Y k,1 Y k,2 |h k ) = 0. (E 4 ) E(Y 2 k,2 |h k ) 4 α 2 L d . (E 5 ) E(Y 4 k,2 |h k ) 2048 s 4 α 4 L d 2 . (E 6 ) E(Y 4 k,1 |h k ) C 3072 s 6 d 2 + 384 s 4 d + 1 α 8 L . (E 7 ) E(Y 2 k,1 |h k ) √ C 128 s 4 d + 1 α 4 L .
For c > 0, we have

P h L 2 h 0 2 c = P ln( h L 2 ) -ln( h 0 2 ) ln(c) = P L-1 k=0 ln(1 + Y k ) ln(c) P L-1 k=0 Y k ln(c) (using ln(1 + x) x for x > -1). Let S = L-1 k=0 Y k -E(Y k |h k ). By (E 1 ) and (E 2 ), L-1 k=0 E(Y k |h k ) Lα 2 L .
So, for c > exp(Lα 2 L ),

P h L 2 h 0 2 c P S ln(c) - L-1 k=0 E(Y k |h k ) P(S ln(c) -Lα 2 L ) P S 2 ln(c) -Lα 2 L 2 E(S 2 ) (ln(c) -Lα 2 L ) 2 (17)
(by Markov's inequality.)

It remains to upper bound E(S 2 ). To this aim, note that

E(S 2 ) = L-1 k=0 E Y k -E(Y k |h k ) 2 L-1 k=0 E Y 2 k 4 Lα 2 L d + 128 √ C Lα 4 L s 4 d + √ CLα 4
L (by (E 3 ), (E 4 ), and (E 7 ))

5 Lα 2 L d .
The last inequality is true for α 2 L 1 √ C(d+128s 4 ) . Therefore, by inequality (17), we obtain, for c > exp(Lα 2 L ),

P h L 2 h 0 2 c 5Lα 2 L d ln(c) -Lα 2 L 2 .
We conclude that, for any δ ∈ (0, 1), with probability at least 1δ,

h L 2 h 0 2 < exp Lα 2 L + 5Lα 2 L dδ .
This shows statement (ii) of the proposition.

Next, to prove statement (i), observe that c > 0,

P h L 2 h 0 2 c = P ln( h L 2 ) -ln( h 0 2 ) ln(c) = P L-1 k=0 ln(1 + Y k ) ln(c) = P L-1 k=0 ln(1 + Y k ) ln(c) and ∀k, Y k - 1 2 + P L-1 k=0 ln(1 + Y k ) ln(c) and ∃k, Y k < - 1 2 .
Using the inequality ln(1 + x) xx 2 for x -1 /2, we obtain

P h L 2 h 0 2 c P L-1 k=0 Y k -Y 2 k ln(c) and ∀k, Y k - 1 2 + P L-1 k=0 ln(1 + Y k ) ln(c) and ∃k, Y k < - 1 2 .
Thus,

P h L 2 h 0 2 c P L-1 k=0 Y k -Y 2 k ln(c) + L-1 k=0 P Y k,2 < - 1 2 . ( 18
)
We handle the two terms above on the right-hand side separately. For the first term, let

Z k = Y k -Y 2 k and S = L-1 k=0 Z k -E(Z k |h k ).
Observe that, by (E 1 )-(E 4 ) and (E 7 ),

L-1 k=0 E(Z k |h k ) Lα 2 L 2 -4 Lα 2 L d -128 √ C Lα 4 L s 4 d - √ CLα 4 L 3 8 Lα 2 L , (19)
where the last inequality is valid for d 64 and α 2 L 1 16 √ C(2s 4 +1) . Hence, for 0 < c < exp( 3Lα 2 L/8),

P L-1 k=0 Y k -Y 2 k ln c = P S ln(c) - L-1 k=0 E(Z k |h k ) P S ln(c) - 3Lα 2 L 8 P S 2 ln(c) - 3Lα 2 L 8 2 E(S 2 ) ln(c) - 3Lα 2 L 8 2
(by Markov's inequality.)

Using the c r -inequality (a + b) n 2 n-1 (a n + b n ) respectively for n = 2 and n = 4, we see that

E(S 2 ) = L-1 k=0 E Z k -E(Z k |h k ) 2 L-1 k=0 E Z 2 k 2 L-1 k=0 E Y 2 k + E Y 4 k 2 L-1 k=0 E(Y 2 k,1 ) + E(Y 2 k,2 ) + 2E(Y k,1 Y k,2 ) + 8E(Y 4 k,1 ) + 8E(Y 4 k,2 ).
By (E 3 )-(E 7 ), it is easy to verify that, for d 64 and

α 2 L 1 ( √ Cs 4 /16+2 √ C+8s 4 )d , E(S 2 ) 10 Lα 2 L d .
This shows that, for c < exp( 3Lα 2 L/8),

P L-1 k=0 Y k -Y 2 k ln(c) 10Lα 2 L d ln(c) - 3Lα 2 L 8 2 .
To conclude the proof, it remains to upper bound the second term of inequality (18). According to inequality (22) in the proof of Lemma 15 (with t = 1 /2), one has

L-1 k=0 P Y k,2 < - 1 2 2L exp - d 64α 2
L s 2 . Putting everything together, we are led to

P h L 2 h 0 2 c 10Lα 2 L d ln(c) - 3Lα 2 L 8 2 + 2L exp - d 64α 2 L s 2 . Take δ ∈ (0, 1). Then, if 2L exp -d 64α 2 L s 2 δ 11 , with probability at least 1 -δ, h L 2 h 0 2 > exp 3Lα 2 L 8 - 11Lα 2 L dδ .
Notice that this inequality is valid under the assumption

α 2 L 2 ( √ Cs 4 +4 √ C+16s 4 )d .

A.4 Proof of Corollary 4

Statement (i) is a consequence of Proposition 2, whereas (ii) is a consequence of Proposition 3 (i). The latter is valid under the conditions d 64 and α L 2 ( √ Cs 4 +4 √ C+16s 4 )d , which is automatically satisfied for all L large enough. Furthermore, an inspection of the proof of Proposition 3 reveals that the divergence in high probability of h L can be proved under the relaxed assumption d 9. Indeed, the main constraint on d comes from the lower bound ( 19), where one needs to make sure that To prove (iii), we use a union bound on both statements of Proposition 3.

A.5 Proof of Proposition 5

The first claim follows from the observation that

∂g(h k , θ k+1 ) ∂h q k =      σ (h k,1 ) 0 . . . 0 0 σ (h k,2 ) . . . 0 . . . . . . . . . . . . 0 0 . . . σ (h k,d )      q k ,
from (A 1 ), and from the assumption on σ .

Let us now prove (ii). In the rest of the proof, the subscript k is ignored to lighten the notation. Observe that

∂g(h, θ) ∂h q = V      σ ( W 1 , h ) 0 . . . 0 0 σ ( W 2 , h ) . . . 0 . . . . . . . . . . . . 0 0 . . . σ ( W d , h )      W q.
Denote by D the matrix in the middle of the right-hand side. Then E ∂g(h, θ) ∂h q 2 h, q = E V DW q 2 |h, q = E DW q 2 |h, q (by (A 1 ))

For model res-2, we have

E ∂g(h, θ) ∂h q 2 h, q = E d i=1 d j=1 W ij q j 2 σ ( W i , h ) h, q .
The conclusion follows from the hypothesis that a σ b and E W q 2 |q = q 2 . For model res-3, we have

E ∂g(h, θ) ∂h q 2 h, q = E d i=1 d j=1 W ij q j 2 1 d j=1 W ij h j 0 h, q .
Since the (W ij ) 1 i,j d are symmetric random variables, we conclude that

E ∂g(h, θ) ∂h q 2 h, q = 1 2 E d i=1 d j=1 W ij q j 2 q = 1 2 E W q 2 |q = q 2 2 .
A.6 Proof of Proposition 6

Letting b = p L / p L , as in Assumption (A 3 ), and taking expectation in (10), we obtain

E p 0 2 p L 2 = E(|b q L (z)| 2 ) = 1 d E( q L (z) 2 ) ( 20
)
(by (A 3 )).

The rest of the proof is similar to the proof of Proposition 2. From (9), we have

q k+1 (z) 2 = q k (z) 2 +α 2 L V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) 2 +2α L q k (z), V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) .
By independence of V k+1 from q k (z) and

∂g(h k ,θ k+1 ) ∂h , E q k (z), V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) = 0.
Next,

E V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) 2 = E E V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) 2 h k , θ k+1 , q k (z) = E ∂g(h k , θ k+1 ) ∂h q k (z) 2 (by (A 1 )) = E E ∂g(h k , θ k+1 ) ∂h q k (z) 2 h k , q k (z) .
By Assumption (A 3 ), we are led to

1 + 1 2 α 2 L E( q k (z) 2 ) E( q k+1 (z) 2 ) (1 + α 2 L )E( q k (z) 2 ),
and thus, by induction, since q 0 (z) = z and E( z 2 ) = d,

d 1 + 1 2 α 2 L k E( q k (z) 2 ) d(1 + 4α 2 L ) k .
In particular, for k = L,

d 1 + 1 2 α 2 L L E( q L (z) 2 ) d(1 + α 2 L ) L .
Therefore, by (20),

1 + 1 2 α 2 L L E p 0 2 p L 2 (1 + α 2 L ) L .
To finish the proof, observe that

1 p L (p 0 -p L ) z = b (q L (z) -z).
Using arguments similar to (20), we may write

E p 0 -p L 2 p L 2 = 1 d E q L (z) -z 2 z 2 . Now, upon noting that q L (z) -z = q L (z) -q 0 (z) = α L L-1 k=0 V k+1 ∂g(h k ,θ k+1 ) ∂h q k (z), E( q L (z) -z 2 ) = α 2 L L-1 k,k =0 E q k (z) ∂g(h k , θ k+1 ) ∂h V k+1 V k +1 ∂g(h k , θ k +1 ) ∂h q k (z) = α 2 L L-1 k=0 E V k+1 ∂g(h k , θ k+1 ) ∂h q k (z) 2 dα 2 L L-1 k=0 (1 + α 2 L ) k = d (1 + α 2 L ) L -1 d exp(Lα 2 L ) -1 2dLα 2 L ,
for Lα 2 L 1. Note that the second equality is obtained by conditioning on every random variable except V k +1 for k < k (and V k+1 for k > k ). Finally, by using Markov's inequality, we conclude that, for any ε > 0,

P p 0 -p L 2 ε p L 2 2Lα 2 L ε .

A.7 Proof of Proposition 7

The proof of Proposition 6 reveals that

E p 0 -p L 2 p L 2 (1 + α 2 L ) L -1.
Using similar arguments, one has

E p 0 -p L 2 p L 2 = 1 d E q L (z) -z 2 z 2 α 2 L L-1 k=0 1 + 1 2 α 2 L k = 1 + 1 2 α 2 L L -1.

A.8 Proof of Corollary 8

The first statement is an immediate consequence of Proposition 6. The second one is a consequence of Proposition 7 and the fact that, for β < 1,

1 + 1 L β L = exp L ln 1 + 1 L β ∼ exp L 1-β → ∞.
Finally, (iii) follows from Proposition 7.

A.9 Proof of Proposition 10

The proposition is a consequence of Kloeden and Platen (1992, Theorems 4.5.3 and 10.2.2) for the SDE

dH t = d 2 σ(H t )dB t .
Letting a(h, t) = 0 and b(h, t) = d 2 σ(h), we need to check the following assumptions:

(H 1 ) The functions a(•, •) and b(•, •) are jointly measurable on R d × [0, 1]. (H 2 ) There exists a constant C 1 > 0 such that, for any x, y ∈ R d , t ∈ [0, 1], a(x, t) -a(y, t) + b(x, t) -b(y, t) C 1 x -y . (H 3 ) There exists a constant C 2 > 0 such that, for any x ∈ R d , t ∈ [0, 1], a(x, t) + b(x, t) C 2 (1 + x ). (H 4 ) E |H 0 2 < ∞. (H 5 ) There exists a constant C 3 > 0 such that, for any x ∈ R d , s, t ∈ [0, 1], a(x, t) -a(x, s) + b(x, t) -b(x, s) C 3 (1 + x )|t -s| 1 /2 .
Assumptions (H 1 ), (H 4 ), and (H 5 ) readily follow from the definitions. Assumption (H 2 ) is true since σ is Lipschitz continuous, and (H 3 ) follows from

σ(x) b x x 1 + x .
A.10 Proof of Proposition 11

Let ψ : R d × [0, 1] → R d be defined for any h ∈ R d , t ∈ [0, 1], by ψ(h, t) = V t g(h, Θ t ).
With this notation, the ODE ( 14) is equivalent to the initial value problem

dH t = ψ(H t , t)dt, H 0 = Ax.
By Assumptions (A 5 ) and (A 6 ), ψ is Lipschitz continuous in its first argument, in the sense that there exists K > 0 such that, for all h

, h ∈ R d , t ∈ [0, 1], ψ(h, t) -ψ(h , t) K h -h .
In addition, it is continuous in its second one. Thus, according to the Picard-Lindelöf theorem (Theorem 22 in Appendix D), this is enough to show that the neural ODE ( 14) has a unique solution on [0, 1]. Note that the solution H is continuous on [0, 1] and is therefore bounded by a constant M > 0.

In order to prove the approximation bound of Proposition 11, we start by proving that both ψ and H are Lipschitz continuous in t. Under (A 5 ) and (A 6 ), this is clear for ψ since H is bounded. Moreover, for any [s, t] ⊂ [0, 1], we have

H t -H s = t s ψ(H u , u)du t s ψ(H u , u) du (t -s) sup u∈[0,1] h∈R d , h M ψ(h, u) .
Now, let K 1 and K 2 denote the Lipschitz constants of ψ (in both arguments) and H respectively, and, for any 0 k L, let t k = k /L. Then we have, for k 1,

H t k -h k = H t k-1 + t k t k-1 ψ(H u , u)du -h k-1 - 1 L ψ(h k-1 , t k-1 ) H t k-1 -h k-1 + t k t k-1 ψ(H u , u) -ψ(h k-1 , t k-1 ) du H t k-1 -h k-1 + K 1 t k t k-1 H u -h k-1 du + K 1 t k t k-1 |u -t k-1 |du 1 + K 1 L H t k-1 -h k-1 + K 1 t k t k-1 H u -H t k-1 du + K 1 t k t k-1 |u -t k-1 |du 1 + K 1 L H t k-1 -h k-1 + (K 2 + 1)K 1 t k t k-1 |u -t k-1 |du = 1 + K 1 L H t k-1 -h k-1 + (K 2 + 1)K 1 2L 2 .
By recurrence, we obtain

H t k -h k k-1 j=0 1 + K 1 L j × (K 2 + 1)K 1 2L 2 L 1 + K 1 L L × (K 2 + 1)K 1 2L 2 e K 1 (K 2 + 1)K 1 2L ,
which concludes the proof.

A.11 Proof of Proposition 12

Starting from (4) and using Assumption (A 6 ), one easily obtains the existence of C 1 and C 2 (whose values depend on the realization of V and Θ) such that

h k+1 (1 + C 1 α L ) h k + C 2 α L .
By recurrence,

h k+1 (1 + C 1 α L ) k h 0 + C 2 C 1 . Hence, using α L 1 /L, h k+1 exp(C 1 ) h 0 + C 2 C 1 .
Since g is Lipschitz continuous on compact sets, it is bounded on every ball of R d × R p . The result is then a consequence of the identity

h L -h 0 = α L L-1 k=0 V k+1 g(h k , θ k+1 ),
since we showed that each term in the sum is bounded by some constant C 3 > 0, independent of L and k. Hence we have that

h L -h 0 C 3 Lα L = C 3 L 1-β ,
yielding the results depending on the value of β.

A.12 Proof of Proposition 13

In the linear case, (4) can be written

h k+1 = h k + α L V k+1 h k , 0 k L -1.
Take y a unit-norm eigenvector of V 0 with associated eigenvalue λ > 0. Then

h k+1 , y = h k + α L V k+1 h k , y = h k , y + α L h k , V k+1 y = h k , y + λα L h k , y + α L h k , (V k+1 -V 0 ) y . Since V is Lipschitz and V k+1 = V k+1/L , there exists c such that V k+1 -V 0 c k+1 L . Hence | h k+1 , y | (1 + λα L )| h k , y | -cα L k + 1 L h k .
Then, by recurrence,

| h L , y | (1 + λα L ) L | h 0 , y | -c α L L L-1 k=0 (k + 1)(1 + λα L ) k h k (1 + λα L ) L | h 0 , y | -cα L (1 + λα L ) L max k h k . Let M = | h 0 ,y | 2cα L , and suppose that h k M for all 0 k L. Then h L | h L , y | (by the Cauchy-Schwartz inequality) (1 + λα L ) L | h 0 , y | -cM α L .
Then, for λα L 1,

h L 1 2 (1 + λα L ) L | h 0 , y | 1 2 exp λLα L 2 | h 0 , y |.
Thus, since Lα L = L 1-β , we have that h L → ∞, which contradicts our assumption that h k M for all 0 k L. We deduce that, for all L large enough,

max k h k > | h 0 , y | 2cα L L→∞ ----→ ∞. Furthermore, max k h k -h 0 h 0 > | h 0 , y | 2c h 0 α L -1 L→∞ ----→ ∞.

B. Technical results

B.1 Lemmas 16 and 17

Lemma 16 Let W ∈ R d×d be a matrix whose entries are symmetric i.i.d. random variables, finite variance, and let σ be an activation function such that, for all

x ∈ R, a|x| |σ(x)| b|x|, 1 / √ 2 a < b 1. Then, for any x ∈ R d , 1 2 E W x 2 E σ(W x) 2 E W x 2 and E ReLU(W x) 2 = 1 2 E W x 2 .
Proof The first part is a consequence of the assumption on σ. To prove the equality, let

X i = d j=1 W ij x j . Then E ReLU(W x) 2 = E d i=1 d j=1 W ij x j 2 1 d j=1 W ij x j 0 = E d i=1 X 2 i 1 X i 0 .
Since the (W ij ) 1 j d are symmetric and independent random variables, X i is also symmetric. Hence E(X 2 i 1 X i 0 ) = 1/2E(X 2 i ), which concludes the proof.

Lemma 17 Let W ∈ R d×d be a matrix whose entries are centered i.i.d. random variables, with finite variance s 2 . Then, for any x ∈ R d , E W x 2 = s 2 d x 2 .

Proof For any 1 i d,

|W x| 2 i = d j=1 W ij x j 2 = d j,j =1 W ij W ij x j x j .
Thus, by independence,

E |W x| 2 i = E d j,j =1 W ij W ij x j x j = d j=1 E(W 2 ij )x 2 j = s 2 x 2 . ( 21
)
The result follows by summing over all i ∈ {1, . . . , d}.

B.2 Proof of Lemma 15

(E 1 ) and (E 2 ) are simple consequences of Assumptions (A 1 ) and (A 2 ).

To prove (E 3 ), let f (h k , θ k+1 ) = V k+1 g(h k , θ k+1 ). Then

E(Y k,2 Y k,1 |h k ) = 1 h k 4 E f (h k , θ k+1 ) 2 h k , f (h k , θ k+1 ) h k = E d i=1 d j=1 f (h k , θ k+1 ) 2 i (h k ) j f (h k , θ k+1 ) j h k .
It is easy to verify that, under Assumption (A 1 ), each term of the sum above has zero expectation. This shows (E 3 ).

To establish (E 4 ), we start by noting that

E h k h k , f (h k , θ k+1 ) h k 2 h k = 1 h k 4 E h k f (h k , θ k+1 )f (h k , θ k+1 ) h k |h k = 1 h k 4 h k E f (h k , θ k+1 )f (h k , θ k+1 ) |h k h k .
Clearly, E(f (h k , θ k+1 ) i f (h k , θ k+1 ) j ) = 0 for i = j. Since, furthermore, the coordinates of f (h k , θ k+1 ) are identically distributed conditionally on h k , we obtain

E f (h k , θ k+1 )f (h k , θ k+1 ) |h k = 1 d E( f (h k , θ k+1 ) 2 |h k )I d .
Thus,

E h k h k , f (h k , θ k+1 ) h k 2 h k = 1 d h k 4 E( f (h k , θ k+1 ) 2 |h k )h k h k 1 d ,
by Assumptions (A 1 ) and (A 2 ). To prove (E 5 ), let ϕ = V k+1 g(h k ,θ k+1 ),h k g(h k ,θ k+1 ) h k . Then, for any t > 0,

P(|Y k,2 | > t) = P |ϕ| > t h k 2α L g(h k , θ k+1 ) = E P |ϕ| > t h k 2α L g(h k , θ k+1 ) h k , θ k+1 .
So, by Lemma 20 in Appendix C,

P(|Y k,2 | > t) E 2 exp - dt 2 h k 2 16α 2 L s 2 g(h k , θ k+1 ) 2 = E E 2 exp - dt 2 h k 2 16α 2 L s 2 g(h k , θ k+1 ) 2 h k E 2 exp - dt 2 h k 2 16α 2 L s 2 E( g(h k , θ k+1 ) 2 |h k ) ,
Proof For any 1 i, j d, set X ij =

x i V ij y j

x y . Let X = R d 2 endowed with the 1 norm, let X be the vector in X whose (id + j)-th coordinate is X ij , and let the function ϕ be defined by

ϕ : X Y -→ d 2 i=1 Y i .
By the triangle inequality, ϕ is a Lipschitz continuous function, with Lipschitz constant equal to 1. Observe also that X ij is a x 2 i s 2 y 2 j/d x 2 y 2 sub-Gaussian. Thus, according to Lemma 19, the sub-Gaussian diameter of X ij is less than Proof The proof is similar to the one of Lemma 20, with X ij =

√ 2x i sy j/ √ d x
V ij x j

x , X = R d , and

ϕ i : X X → d j=1 X ij .
Each function ϕ i is a Lipschitz continuous function, with Lipschitz constant equal to 1. Observe now that the random variable X ij is x 2 j s 2 /d x 2 sub-Gaussian. Thus, according to Lemma 19, the sub-Gaussian diameter of X ij is less than √ 2x j s / √ d x . Therefore, according to Kontorovich (2014, Theorem 1), for any t > 0, P (ϕ i (X) t) 2 exp -

t 2 2 d j=1 2s 2 x 2 j d x 2
, that is

P | V i , x | x t 2 exp - dt 2 4s 2 .
Hence (see, e.g., [START_REF] Pauwels | Statistics, Optimization and Algorithms in High Dimension[END_REF],

E V i , x x 2q q! 8s 2 d q . ( 23
)
From identity (21) in the proof of technical Lemma 17, given in Appendix B, we obtain that, for q = 1,

E V i , x x 2 = 1 d , (24)
which is an improvement by a factor 8s 2 over the previous upper bound. To conclude, it remains to conclude V x 4 and V x 8 with the V i , x . To do so, observe that

V x 4 = d i=1 V i , x 2 2 = d i,j=1 i =j V i , x 2 V j , x 2 + d i=1 V i , x 4 .
Hence, by independence of the (V i ) 23) and ( 24)) Similarly,

1 i d , E V x 4 x 4 = d i,j=1 i =j E V i , x 2 x 2 E V j , x 2 x 2 + d i=1 E V i , x 4 x 4 = d(d -1) 1 d 2 + d 2(8s 2 ) 2 d 2 1 + 128s 4 d (by (
V x 8 = d i=1 V i , x 2 3 = d i,j,k=1 i =j =k V i , x 2 V j , x 2 V j , x 2 + d i,j=1 i =j V i , x 2 V j , x 4 + d i=1 V i , x 8 .
Hence,

E V x 8 x 8 = d i,j,k=1 i =j =k E V i , x 2 x 2 E V j , x 2 x 2 E V k , x 2 x 2 + d i,j=1 i =j E V i , x 4 x 4 E V j , x 2 x 2 + d i=1 E V i ,
Ψ(Y ) v -Ψ(Y ) v v s f (Y u , u) -f (Y u , u) du v s K f Y u -Y u du K f v s Y -Y ∞ du K f Y -Y ∞ (t -s).
This yields

Ψ(Y ) -Ψ(Y ) ∞ K f (t -s) Y -Y ∞ ,
which means that the function Ψ is Lipschitz continuous on C ([s, t], R d ) endowed with the supremum norm, with Lipschitz constant K f (ts). So, on any interval [s, t] of length smaller than δ = 1 /2K f , the function Ψ is a contraction. Thus, by the Banach fixed-point theorem, for any initial value ζ, Ψ has a unique fixed point. Hence, there exists a unique solution to (25) on any interval of length δ with any initial condition. To obtain a solution on [0, 1] it is sufficient to concatenate these solutions.

E. Detailed experimental setting

Our code is available at https://github.com/PierreMarion23/scaling-resnets.

To obtain Figures 1 to 3, we initialize ResNets from res-3 with the hyper-parameters of Table 2.

Each experiment is repeated 50 times, with independent data and weight sampling. For Figures 4 and5, we take the same hyper-parameters except for β, which now takes values in {0.5, 1, 2}, and for the weight distribution. The weights are now initialized as

Figure 1 :

 1 Figure 1: Evolution of h Lh 0 / h 0 as a function of L for different values of β and an i.i.d. U(-3/d, 3/d) initialization of model res-3, with d = 40. The input is a random Gaussian observation x in dimension n in = 64. The experiment is repeated with 50 independent randomizations.

Figure 2 :

 2 Figure 2: Empirical distributions of the norms for β = 1 /2, L = 10 3 , d = 100. The experiment is repeated with 10 4 independent randomizations.

Figure 3 :

 3 Figure 3: Evolution of p 0p L / p L as a function of L for different values of β and an i.i.d. U(-3/d, 3/d) initialization of model res-3, with d = 40. The input is a random Gaussian observation x in dimension n in = 64. The experiment is repeated with 50 independent randomizations.

Figure 4 :

 4 Figure 4: Evolution of h Lh 0 / h 0 as a function of L for different values of β and a smooth initialization of model res-3, with d = 40. The input is a random Gaussian observation x in dimension n in = 64. The experiment is repeated with 50 independent randomizations.

Figure 5 :

 5 Figure 5: Evolution of p 0p L / p L as a function of L for different values of β and a smooth initialization of model res-3, with d = 40. The input is a random Gaussian observation x in dimension n in = 64. The experiment is repeated with 50 independent randomizations.

Figure 6 :

 6 Figure 6: Examples of realizations of a fractional Brownian motion B H for different Hurst indexes H. Note that the smaller the value of H, the more irregular the trajectory is.

  ln 10 ( p 0p L / p L )

  Figure 8: Plot of a given coordinate of θ k , after training, as a function of the layer index k ranging from 1 to the depth L = 1000 for three different choices of β and initializations.

Figure 9 :

 9 Figure 9: Accuracy after training as a function of the regularity of the weights at initialization and scaling. For each point of the heatmap, the model was trained on a grid of learning rates, and the best performance is shown.

  which is the case for d = 9.

  y . By Kontorovich (2014, Theorem 1), for any t > 0, one hasP (ϕ(X) t) 2 exp -Lemma 21 (Bound of moments of quadratic forms) Let V be a R d×d matrix whose entries are i.i.d s 2 /d sub-Gaussian random variables, with variance 1 /d. Then, for any x ∈ R d ,

						t 2 i,j=1 d x 2 y 2 2 d 2s 2 x 2 j i y 2	,
	that is		P	V x, y x y	t	2 exp -	dt 2 4s 2 .
	x = 0,
	E	V x 4 x 4	1 +	128s 4 d	and E	V x 8 x 8	1 +	384s 4 d	+	3072s 6 d 2 .

  . A version of the Picard-Lindelöf theorem Theorem 22 Assume that f is Lipschitz continuous in its first argument and continuous in its second one. Then, for any z ∈ R d , the initial value problemdH t = f (H t , t)dt, H 0 = z, (25)admits a unique solution H : [0, 1] → R d . Proof Let C ([s, t], R d ) be the set of continuous functions from [s, t] to R d . For any [s, t] ⊂ [0, 1], ζ ∈ R d , let Ψ be the function Ψ : C ([s, t], R d ) → C ([s, t], R d ) For any Y, Y ∈ C ([s, t], R d ), v ∈ [s, t], one has, denoting by K f the Lipschitz constant of f in its first argument,

					x 8
					x 8
	= d(d -1)(d -2)	1 d 2 + 3d(d -1)	2(8s 2 ) 2 d 3 + d	6(8s 2 ) 3 d 3
	1 +	384s 4 d	+	3072s 6 d 2 .

DY → v → ζ + v s f (Y u , u)du .

Acknowledgments

The authors thank S. Schoenholz for fruitful discussion. P. Marion has been supported by a grant from Région Île-de-France and by a Google PhD Fellowship award.

A. Proofs

Throughout the proofs, the i-th coordinate of a vector v is denoted by v i . Similarly, the i-th row of a matrix M is denoted by M i , and its (i, j)-th entry by M ij .

A.1 Proof of Proposition 1

Statement (i) is clear (with C = 1) since, for any h ∈ R d ,

by Jensen's inequality. Finally, using Assumption (A 2 ), we deduce that

In particular, for all q 1 (see, e.g., [START_REF] Pauwels | Statistics, Optimization and Algorithms in High Dimension[END_REF],

The result is obtained by taking q = 2. Finally, (E 6 ) and (E 7 ) are consequences of Lemma 21 in Appendix C.

C. Concentration of sub-Gaussian random matrices

In this appendix, we are interested in the concentration of linear and quadratic forms of sub-Gaussian matrices (Lemma 20 and Lemma 21). These two propositions are byproducts of the main result of [START_REF] Kontorovich | Concentration in unbounded metric spaces and algorithmic stability[END_REF], which generalizes McDiarmid's inequality to sub-Gaussian variables. We start by a technical result regarding the sub-Gaussian diameter introduced by Kontorovich ( 2014), whose definition is recalled below.

Definition 18 Let X be a real-valued random variable, X an independent copy of X, and ε a Rademacher random variable, independent of X and X . Then the sub-Gaussian diameter of X is defined as the smallest t such that ε|X -X | is t 2 sub-Gaussian.

Lemma 19 Let X be a s 2 sub-Gaussian symmetric random variable. Then the sub-Gaussian diameter of X is less than √ 2s.

Proof Let λ ∈ R. Then, using the notation of Definition 18, one has

where the last equality is a consequence of the symmetry of X.

We are now ready to prove the two main results of this appendix.

Lemma 20 (Bound on the deviation of linear forms) Let V be a R d×d matrix whose entries are i.i.d s 2 /d sub-Gaussian random variables. Then, for any x, y ∈ R d , x, y = 0, 7 discretizations of a Gaussian process. More precisely, each entry of V and Θ is an independent Gaussian process with zero mean and an RBF kernel of variance 10 -2 .

To obtain Figure 7, we take the hyper-parameters of Table 3.

More precisely, for each 1 i, j d, we let (V k+1,i,j ) 0 k L-1 be the increments of a fractional Brownian motion (fBm), where the various fBm involved are independent. The procedure is the same for θ.

In Figure 9, we use res-1, with the hyper-parameters of We train on MNIST 1 and CIFAR-10 2 using the Adam optimizer [START_REF] Kingma | A method for stochastic optimization[END_REF]) for 10 epochs. The learning rate is divided by 10 after 5 epochs. The best performance on the learning rate grid is reported in the figure .  Figure 8 is obtained by plotting a random coordinate of θ k , after training on MNIST.