Carlos Esteve-Yagüe

Borjan Geshkovski

AND Dario Pighin

Enrique Zuazua

Large-time asymptotics in deep learning

Keywords: Deep learning, ResNets, neural ODEs, regularization path, optimal control, exponential stability, approximation, turnpike theory AMS Subject Classification. 49J15, 49M15, 49J20, 49K20, 93C20, 49N05

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Modern supervised learning addresses the problem of predicting from data, which roughly consists in approximating an unknown function f : X → Y from N known but possibly noisy samples { x i , y i } N i=1 ⊂ X ×Y. Depending on the nature of the labels y i , one distinguishes two types of supervised learning tasks, namely that of classification (labels take values in a finite set of m classes, e.g. Y := {1, . . . , m}) and regression (labels take continuous values in Y ⊂ R m ). In many applications, the dimension d of each sample x i ∈ X ⊂ R d may be big compared to the number/dimension m of the labels -in image classification for instance, a sample of the ImageNet dataset [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], which has m = 1000 classes, is an element of R 65536 .

A plethora of methods for finding f (•) efficiently with theoretical and empirical guarantees have been developed and investigated in the machine learning literature in recent decades. Prominent examples, to name a few, include linear parametric methods (e.g. linear or logistic regression), kernel-based methods (e.g. support vector machines), treebased methods (e.g. decision trees) and so on. We refer to [START_REF] Goodfellow | Deep learning[END_REF] for a comprehensive presentation of these topics.

Deep neural networks are parametrized computational architectures which propagate each individual sample of the input data { x i } N i=1 ⊂ X ⊂ R d across a sequence of linear parametric operators and simple nonlinearities. The so-called residual neural networks (ResNets, [START_REF] He | Deep residual learning for image recognition[END_REF]) may, in the simplest case, be cast as schemes of the mould

   x k+1 i = x k i + σ w k x k i + b k for k ∈ {0, . . . , N layers -1} x 0 i = x i ∈ R d (1.1)
for all i ∈ [N ], where we set [N ] := {1, . . . , N }. The unknowns are the states x k i ∈ R d for any i ∈ [N ], σ is an explicit scalar, Lipschitz continuous nonlinear function defined component-wise in (1.1), w k , b k N layers -1 k=0 are optimizable parameters (controls) with w k ∈ R d×d -called weights, and b k ∈ R d -called biases, and N layers 1 designates the number of layers referred to as the depth. The training process consists in finding optimal parameters steering all of the network outputs x N layers i as close as possible to the corresponding labels y i by solving

min {w k ,b k } N layers -1 k=0 1 N N i=1 loss P x N layers i
, y i , whilst guaranteeing reliable performance on unseen data (ensuring generalization). Here loss(•, •) is a given continuous and nonnegative function which differs depending on the task in hand -for instance loss(x, y) := xy p p for p ∈ {1, 2} is commonly used for regression tasks, while loss(x, y) = log(1 + exp(-yx)) may be used for binary classification, namely when y i ∈ {-1, 1} (we refer to Section 2.3 for more general settings). On the other hand, P : R d → R m is an affine map whose coefficients in practice are part of the optimizable parameters. In our work, we shall assume that P is given and specified on a case-by-case basis.

Due to the inherent dynamical systems nature of ResNets, several recent works have aimed at studying an associated continuous-time formulation in some detail, a trend started with the works [START_REF] Weinan | A proposal on machine learning via dynamical systems[END_REF][START_REF] Haber | Stable architectures for deep neural networks[END_REF] (see also [START_REF] Lu | Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations[END_REF][START_REF] Ruthotto | Deep neural networks motivated by partial differential equations[END_REF]). This perspective is motivated by the simple observation that for any i ∈ [N ] and for T > 0, (1.1) is roughly the forward Euler scheme for the neural ordinary differential equation (neural ODE) ẋi (t) = σ(w(t)x i (t) + b(t)) for t ∈ (0, T )

x i (0) = x i ∈ R d .
(1.2)

The continuous-time, neural ODE formalism of deep learning has been used to great effect in applications -for instance, by using adaptive ODE solvers [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Dupont | Augmented Neural ODEs[END_REF][START_REF] Queiruga | Continuousin-depth neural networks[END_REF], symplectic and multigrid methods [START_REF] Celledoni | Structure preserving deep learning[END_REF][START_REF] Gunther | Layer-parallel training of deep residual neural networks[END_REF], or indirect training algorithms based on the Pontryagin Maximum Principle [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF][START_REF] Benning | Deep learning as optimal control problems: Models and numerical methods[END_REF] -, and also for generative modeling through normalizing flows [START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF][START_REF] Chen | Residual flows for invertible generative modeling[END_REF]. We emphasize that the origins of continuous-time supervised learning go back to the 1980s -the neural network model proposed in [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] is a differential equation, whereas in [START_REF] Lecun | A theoretical framework for back-propagation[END_REF] back-propagation is connected to the adjoint method arising in optimal control. Related works include studies on identification of the weights from data [Albertini and Sontag, 1993;Albertini et al., 1993] and controllability of continuous-time recurrent networks [START_REF] Sontag | Complete controllability of continuous-time recurrent neural networks[END_REF][START_REF] Sontag | Further results on controllability of recurrent neural networks[END_REF]. The role of the final time horizon T > 0, which plays a key role in the control of dynamical systems, has not been analyzed in the context of supervised learning problems via models such as (1.2). Now note once again that in the ResNet (1.1) the time-step ∆t = T /Nlayers is fixed (equal to 1), and each time instance of a forward Euler discretization to (1.2) would represent a different layer of (1.1). Hence, whenever the time-step ∆t = T /Nlayers is fixed (or goes to zero when T is increased), the time horizon T in (1.2) serves as an indicator of the number of layers N layers in the ResNet (1.1). Thus, a good knowledge of the behavior of the learning problem and the neural ODE flow over longer time horizons is desirable in view of understanding approximation and generalization properties. In this work, we aim to bridge this gap by proposing several insights stemming from an analysis of the role of the time horizon T .

1.1. Our contributions. We shall focus this presentation on the neural ODE (1.2), but our results also hold for other systems, as seen in the respective statements.

1. We first consider the classical regularized empirical risk minimization problem inf Let us put the above result into context. For neural ODEs for which L 2regularization suffices, we remark that T -→ ∞ is equivalent to λ 0. The latter is the regularization path limit, studied in the literature for linear models and multi-layer perceptrons (but not for more compound neural ODE models), where the solutions obtained in the limit can be shown to satisfy desirable generalization properties (see Section 1.2).

Using similar arguments as when T -→ ∞, in Theorem 3.2 we obtain the same conclusions when λ 0 and T is fixed. Consequently, Theorem 3.1 stipulates generalization properties -namely, optimizing with T 1, which may be interpreted as a larger depth for ResNets, has the practically desirable effect of making the training error close to zero, but by means of parameters with the smallest amplitude and thus trajectories with the least oscillations. 2. In the setting of losses for which E does not attain its minimum equal to zero, occurring in many classification contexts (for instance loss(P x, y) = log 1 + e -yP x when y i ∈ Y = {-1, 1} or multi-label tasks where y i ∈ [m] for m 2 via crossentropy loss), we show that the training error E(x T (T )) is at most of the order O 1 T α for all α ∈ (0, 1) (see Theorem 3.3). Analog results hold when T > 0 is fixed and λ 0 (see Theorem 3.4). 3. To enhance the polynomial stability of the training error E(x T (T )) with respect to T , we also consider the augmented empirical risk minimization problem 2 inf [w,b]∈H(0,T ;R du )

x i (•) solves (1.2) E(x(T )) + 1 N N i=1 T 0 x i (t) -x i 2 dt + λ [w, b] 2 H(0,T ;R du ) , (1.5)
where P : R d → R m appearing in E is Lipschitz and surjective, loss is an pdistance, while the targets x i ∈ P -1 ({ y i }) for all i ∈ [N ] are arbitrary but given. Under a specific controllability assumption but without any differentiability assumptions on the activation function σ or smallness assumptions on the dataset, in Theorem 4.1 (see also Theorem 4.2) we show that optimal parameters [w T (t), b T (t)] and the training error E(x T (t)) of the corresponding vector x T (t) of solutions to (1.1) are at most of the order O(e -µt ) for any t ∈ [0, T ] and some µ > 0 independent of T . This result is in line with Theorem 3.1 and Theorem 3.3, but with an improved estimate of the time horizon needed to be ε-close to the interpolation regime for any given ε > 0. Due to the exponentially small global minimizers, numerical experiments show that the learned flow is simple, which could stipulate possible generalization properties. Theorem 4.1 is a manifestation of the 2 As for (1.3), we consider H = L 2 when existence of minimizers can be ensured; however, unlike for (1.3), in the setting of (1.2), we will sometimes consider a BV regularization (see Section 1.4 for a definition) rather than the stronger H 1 regularization for technical reasons (see Remark 6). so-called turnpike property 3 , well-known in optimal control theory ( [START_REF] Porretta | Long time versus steady state optimal control[END_REF][START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF][START_REF] Faulwasser | Turnpike properties in optimal control: An overview of discrete-time and continuous-time results[END_REF]).

Problem (1.5) is motivated by the natural training problem which aims at minimizing the error over every time/layer:

inf [w,b]∈H(0,T ;R du ) x i (•) solves (1.2) T 0 E(x(t)) dt + λ [w, b] 2 H(0,T ;R du ) ,
(1.6)

where loss(•, •) appearing in E is continuous and nonnegative, but otherwise arbitrary. Whilst left without proof, numerical experiments stipulate that a similar decay for the training error, and, combined with Theorem 4.1, motivate the usage of (1.6) in practice (see Section 4.1). 4. In Theorem 5.1 we show that the controllability assumption needed for Theorem 4.1 is satisfied for a subclass of neural ODEs with C 1 -regular activation functions σ. We also illustrate the relation between the amplitude of the weights needed to reach the interpolation regime, and the dispersion of the input data (see Proposition 5.1). 5. To address variable width architectures motivated by multi-layer perceptrons and convolutional neural networks, in Section 6 we study a continuous spacetime neural network formulation taking the form of a scalar integro-differential equation proposed in [START_REF] Liu | Selection dynamics for deep neural networks[END_REF]. We show that, by making use of a time-dependent moving grid for discretizing the spatial variable, these equations provide a framework for addressing ResNets with variable widths. Furthermore, in Theorem 6.1 (resp. Theorem 6.2), we show that some of our finite-dimensional conclusions from Theorem 3.1 (resp. Theorem 4.1) transfer, under similar assumptions, to the continuous space-time neural networks.

1.2. Related work. We note that the impact of the time horizon for the regularized empirical risk minimization problem has been studied from a computational perspective in [START_REF] Effland | Variational networks: An optimal control approach to early stopping variational methods for image restoration[END_REF]. We also refer the reader to the recent work [START_REF] Faulwasser | On the turnpike to design of deep neural nets: Explicit depth bounds[END_REF] for a stability analysis for the augmented supervised learning problem (1.5) in the ResNet setting, and to [START_REF] Yagüe | Sparse approximation in learning via neural ODEs[END_REF] for a rigorous analysis of the problem (1.6) with an L 1 -regularization of the parameters, where sparsity patterns of the parameters are shown.

In [START_REF] Thorpe | Deep limits of residual neural networks[END_REF] (see also [START_REF] Avelin | Neural ODEs as the deep limit of ResNets with constant weights[END_REF]), the authors show, via Γ-convergence arguments, that the optimal control parameters in the discretetime context converge to those of the continuous-time context when the step-size ∆t goes to 0. The latter is interpreted as an infinite layer limit when the final time horizon T in the continuous-time context is fixed (equal to 1). Our results are of different nature. Rather than proving that the discrete-time controls converge to the continuous-time ones, we exhibit the continuous-time, neural ODE representation, for which the final time horizon indicates the number of layers for the associated time-discretization when the time-step ∆t is fixed, and aim to characterize the impact of T on the optimal parameters and on the training error.

3 The turnpike property is a staple of optimal control problems such as the linear quadratic regulator (LQR) over large but finite time horizons; it refers to the fact that optimal controls (parameters) and trajectories, in long time intervals, remain close to the optimal stationary controls and trajectories. In the setting of (1.5), due to the structure of the neural ODE, the optimal stationary parameters are zero, and thus the optimal stationary state is the collection {xi} i∈[N ] , as reflected in our result.

Our results are also related to several questions studied in existing literature.

Universal approximation. On a first note, the asymptotic results presented herein may (heuristically) be interpreted as approximation results in the sense of the universal approximation theory. These are density results for neural networks, and in the simplest cases can be interpreted in terms of the elementary building blocks of measure theory such as the density of simple functions in Lebesgue spaces. The first result in this direction is the seminal work [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF], which indicates that shallow neural networks with increasing width, i.e. a superposition of sufficiently many dilated and translated sigmoids, may approximate any continuous function on compact sets. We also refer to [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF][START_REF] Pinkus | Approximation theory of the MLP model in neural networks[END_REF] for an extension to multi-layer neural networks. Our results are somewhat dual to [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF] -therein, to increase the approximation accuracy, the width is allowed to grow, whilst we fix the width and allow the depth to increase. We do note however that we prove approximation properties for the optimal parameters, and for a fixed dataset, unlike what is commonly done in universal approximation theorem, where the parameters are not known explicitly. We refer to the thesis [START_REF] Müller | Universal flow approximation with deep residual networks[END_REF] for results and a comprehensive review of universal approximation results for ResNets, and to the recent works [START_REF] Li | Deep learning via dynamical systems: An approximation perspective[END_REF][START_REF] Teshima | Universal approximation property of neural ordinary differential equations[END_REF][START_REF] Zhang | Approximation capabilities of neural ODEs and invertible residual networks[END_REF] and [START_REF] Ruiz-Balet | Neural ODE control for classification, approximation and transport[END_REF], for universal approximation results for neural ODE and for an analysis on the latter's working mechanisms.

Regularization path limit: λ 0. The regularization path limit λ 0 for minimizers [w λ , b λ ] has also been addressed in the machine learning literature, albeit for more classical models. For instance in [START_REF] Rosset | Boosting as a regularized path to a maximum margin classifier[END_REF][START_REF] Rosset | Margin maximizing loss functions[END_REF], the authors study linear logistic regression, and show convergence of the margin to the max-margin as λ 0, assuming linearly separable data. The max-margin, support vector machine solution (see Remark 4 for details) is a special example among all solutions that fit the training data -another example includes minimal 2 -norm solution for linear regression (or generally supervised learning tasks in which the loss may attain its minimum). Both these solutions can be shown to ensure generalization by virtue of explicit generalization error estimates [START_REF] Bartlett | Rademacher and Gaussian complexities: Risk bounds and structural results[END_REF][START_REF] Kakade | On the complexity of linear prediction: Risk bounds, margin bounds, and regularization[END_REF] for linear models or multi-layer perceptrons. This insight could stipulate a likely generalization capacity of our asymptotic limits as T -→ ∞ or λ 0. The results of [START_REF] Rosset | Boosting as a regularized path to a maximum margin classifier[END_REF][START_REF] Rosset | Margin maximizing loss functions[END_REF] have subsequently been extended in [START_REF] Wei | Regularization matters: Generalization and optimization of neural nets vs their induced kernel[END_REF] (and some of the references therein) to perceptrons with ReLU activations, where the intrinsic homogeneity of the network is used -the authors prove an analog result to Theorem 3.2 when λ 0 for regression tasks (the loss is a distance) and two-layer perceptrons. Our results can be seen as an extension of the aforementioned insights to the neural ODE context.

Implicit regularization of gradient methods. When training without explicit regularization (i.e. λ = 0), a common approach in the literature is to resort to algorithmdependent generalization analysis, where the limit solutions are akin to the limit solutions when λ 0. Whilst the former is not a question we address in this work, where generally estimates are provided in terms of the number of algorithm iterations rather than the number of layers, we provide a brief overview of the literature for completeness. We refer the reader to the recent review [START_REF] Bartlett | Deep learning: a statistical viewpoint[END_REF] for a comprehensive presentation of these artifacts.

The implicit bias of gradient descent [START_REF] Soudry | The implicit bias of gradient descent on separable data[END_REF][START_REF] Gunasekar | Implicit bias of gradient descent on linear convolutional networks[END_REF]) indicates that in the overparametrized regime, after training a linear model or perceptron with gradient-based methods until zero training error, without requiring any explicit parameter regularization, among the many predictors which overfit on the training dataset, the algorithm selects the one which performs best on the test dataset (e.g. minimal 2 -norm solution or max-margin solution). Thus, even-though the approach is significantly different from that when λ 0, the asymptotic limits oftentimes coincide. Recent works have shown that gradient descent can allow overparametrized multilayer networks to attain arbitrarily low training error on fairly general datasets ( [START_REF] Du | Gradient descent finds global minima of deep neural networks[END_REF]Allen-Zhu et al., 2019a,b]), and find minimum-norm/maximum-margin solutions that fit the data in the settings of logistic regression, deep linear networks, and symmetric matrix factorization [START_REF] Gunasekar | Implicit bias of gradient descent on linear convolutional networks[END_REF][START_REF] Soudry | The implicit bias of gradient descent on separable data[END_REF][START_REF] Ji | Risk and parameter convergence of logistic regression[END_REF][START_REF] Ma | Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion, and blind deconvolution[END_REF]). In [Chizat andBach, 2018, 2020] overparametrization is approached from the point of view of the width of the neural network, unlike our depth-inspired perspective. The authors consider a two-layer perceptron with ReLU activation, exhibit a mean-field, Wasserstein gradient flow formulation of the descent scheme yielding the optimal parameters, and prove that these parameters approach global minimizers of the cost functional, with the global minimizer being characterized as a max-margin classifier in a certain non-Hilbert space of functions. We refer the reader to [START_REF] Sirignano | Mean field analysis of neural networks: A law of large numbers[END_REF][START_REF] Javanmard | Analysis of a two-layer neural network via displacement convexity[END_REF][START_REF] Mei | A mean field view of the landscape of two-layer neural networks[END_REF][START_REF] Nitanda | Stochastic particle gradient descent for infinite ensembles[END_REF]Chizat et al., 2018] for related works in this direction, and [START_REF] Nguyen | A rigorous framework for the mean field limit of multilayer neural networks[END_REF] for an extension of the aforementioned convergence results to multi-layer perceptrons.

1.3. Outline. The remainder of the paper is organized as follows. In Section 2, we give a comprehensive presentation of the neural ODE perspective of deep supervised learning. In Section 3, we present our main results in the context of regularized empirical risk minimization (Theorem 3.1 and Theorem 3.3). In Section 4, we present our main result for the augmented empirical risk minimization problem, namely exponential decay of the training error with exponentially small parameters (Theorem 4.1 and also Theorem 4.2), as well as several numerical experiments stimulating the validity of our conjectures. Finally in Section 6, we present the continuous space-time analog of residual neural networks with variable widths, depict some possible approaches for passing from the continuous to the discrete case, and present extensions of Theorem 3.1 and Theorem 4.1 in this context. The proofs of all results may be found in Section 7.

1.4. Notation, conventions, assumptions. We set ẋ(t) := dx dt (t) and [n] := {1, . . . , n}. Given a ∈ R n , we denote by a its transpose. We use the notation x T to display the dependence of a vector x ∈ R n on the time horizon T . We denote by a the standard euclidean norm when a ∈ R n is a vector, and the entry-wise euclidean norm (Frobenius norm) when a ∈ R n×m is a matrix -we recall

a :=   n j=1 m k=1 |a j,k | 2   1 /2 .
We denote by Lip(R) the set of functions f : R -→ R which are globally Lipschitz continuous, and by L 2 (0, T ; R n ) (resp. H 1 (0, T ; R n )) the Lebesgue (resp. Sobolev) space consisting of all functions f : (0, T ) -→ R n which are square integrable (resp. square integrable and with a square integrable weak derivative) -recall that

H 1 (0, T ; R n ) is endowed with the norm f 2 H 1 (0,T ;R n ) := f 2 L 2 (0,T ;R n ) + ḟ 2 L 2 (0,T ;R n ) .
We will use the convention H 0 = L 2 . We recall (see [Ambrosio et al., 2000, p. 117]) the definition of the space of bounded variation functions BV(0, T ; R n ) as the space of integrable functions whose weak derivative is a finite Radon measure:

BV(0, T ; R n ) := f ∈ L 1 (0, T ; R n ) : ∃ Df ∈ M(0, T ; R n ) such that n j=1 T 0 f j (t)ϕ j (t) dt = - n j=1 T 0 ϕ j (t) dDf j (t), ∀ϕ ∈ C 1 c (0, T ; R n ) .
Here, M(0, T ; R n ) denotes the set of finite Radon measures on (0, T ) with values in R n .

The space BV(0, T ; R n ) is endowed with the norm

f BV(0,T ;R n ) := f L 1 (0,T ;R n ) + n j=1 |Df j |(0, T ).
These definitions apply to matrix valued functions by simply considering the vectorized form of the matrix. Finally, we say that a function f :

R → R is α-homogeneous, for α > 0, if f (cx) = c α x for c ∈ R and x ∈ R.
Throughout the remainder of this work, we will work under the following couple of assumptions, which are universal in the context of machine learning.

Assumption 1. We henceforth assume that we are given a training dataset

{ x i , y i } i∈[N ] ⊂ X × Y,
where X ⊂ R d and x i = x j for i = j. We suppose that Y = ∅, with Y ⊂ R m for regression tasks, Y = {-1, 1} for binary classification and Y = [m] for mutli-label classification tasks.

We henceforth denote

d u := d × (d + 1), d x := dN.
We will, most importantly, make use of the following abuse of notation: given the input data { x i } i∈[N ] with x i ∈ R d , we define the stacked vector x 0 ∈ R dx by

x 0 :=    x 1 . . . x N    with x 0 i = x i for i ∈ [N ].
(1.7)

We shall apply convention (1.7) to other vectors generally denoted by x ∈ R dx and defined by stacking given sub-vectors

{x i } i∈[N ] with x i ∈ R d : x :=    x 1 . . . x N    .
Thus, when we write x i for i ∈ [N ], we insinuate that x i ∈ R d , rather than considering solely the i-th row of x ∈ R dx (which is a scalar).

The following assumption is satisfied by most of the commonly used activation functions, including sigmoids such as σ(x) = tanh(x), and rectifiers: σ(x) = max{αx, x} for α ∈ [0, 1).

Assumption 2. Unless stated otherwise, we fix an activation function σ satisfying σ ∈ Lip(R) and σ(0) = 0.

Roadmap to learning via neural ODEs

2.1. Feed-forward neural networks. The canonical example of a feed-forward neural network is the multi-layer perceptron (MLP), which generally takes the form

   x k+1 i = σ w k x k + b k for k ∈ {0, . . . , N layers -1} x 0 i = x i ∈ R d (2.1) for i ∈ [N ];
here N layers 1 is the depth of (2.1), each k is a layer, the vector x k i ∈ R d k designates the state at the layer k, d k is referred to as the width of the layer k, while w k ∈ R d k+1 ×d k and b k ∈ R d k are the optimizable weight and bias parameters of (2.1). Finally, σ ∈ Lip(R) is a fixed nonlinear activation function -by abuse of notation, we define the vector-valued analog of σ component-wise, namely, σ : R d -→ R d is defined by σ(x) j := σ(x j ) for j ∈ [d]. Common choices include sigmoids such as σ(x) = tanh(x) or σ(x) = 1 1+e -x , and rectifiers: σ(x) = max{x, ax} for a fixed 0 a < 1. In practice, the activation σ is generally selected using cross-validation. It can readily be seen that the formulation (2.1) coincides with the more conventional formulation of neural networks as compositional structures of parametric affine operators and nonlinearities, as namely x

N layers i = (σ • Λ k • . . . • σ • Λ 0 )( x i ), with Λ k x := w k x + b k for k ∈ {0, . . . , N layers }.
Note that the iterative nature of the MLP (2.1) stimulates permuting the order of the parametric affine maps and the nonlinearity σ, to the effect of considering the almost equivalent, but somewhat simpler system

   x k+1 i = w k σ x k + b k for k ∈ {0, . . . , N layers -1} x 0 i = x i ∈ R d .
(2.2)

We will henceforth concentrate on residual neural networks (ResNets). Contrary to the multi-layer perceptrons (2.1) -(2.2), when considering ResNets one typically needs to assume that the width d k is fixed over every layer k, namely d k = d for every k. We refer to Section 6 for variable width ResNets. In the fixed width context, a residual neural network generally takes the form d+1) and f as the right hand side in (2.1) or (2.2). As explained in [START_REF] Lu | Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations[END_REF], other network architectures (including convolutional neural networks) can be fit into the residual network framework.

   x k+1 i = x k i + f u k , x k i for k ∈ {0, . . . , N layers -1} x 0 i = x i ∈ R d (2.3) for i ∈ [N ], where x k i ∈ R d for any i, k, u k := [w k , b k ] ∈ R d×(
2.2. Neural ODEs. It is readily seen that (2.3) corresponds, modulo a scaling factor ∆t = T /Nlayers = 1, to the forward Euler discretization of ẋi (t) = f(u(t), x i (t)) in (0, T )

x i (0) = x i ∈ R d , (2.4) for i ∈ [N ].
Here T > 0 is a fixed time horizon, and u(t d+1) . As per what precedes, the nonlinearity f in (2.4) may take the form

) := [w(t), b(t)] ∈ R d×(
f(u(t), x i (t)) := σ(w(t)x i (t) + b(t)) (2.5) or f(u(t), x i (t)) = w(t)σ(x i (t)) + b(t).
(2.6) for i ∈ [N ]. We will address both cases in our analytical study, and emphasize the stark differences between the two. The above parametrizations are not the lone considered in practice. In fact, one may consider, for instance, combinations of (2.5) and (2.6) which allow intermediate exploration (bottlenecks) in different dimensions:

f(u(t), x i (t)) := w 2 (t)σ(w 1 (t)x i (t) + b 1 (t)) + b 2 (t)
(2.7)

where now w

1 (t) ∈ R d hid ×d , w 2 (t) ∈ R d×d hid , b 1 (t) ∈ R d hid and b 2 (t) ∈ R d .
2.3. Training. For an input sample x i ∈ R d , the output of the neural ODE (2.4), which is used for comparison with the corresponding label y i , is a projection of the form P x i (T ) ∈ R m for some affine map P : R d → R m :

P x := p 1 x + p 2 ,
where p 1 ∈ R m×d and p 2 ∈ R m . In the context of binary classification, namely m = 1 with y i = ±1, one may also use P x := tanh(p 1 x + p 2 ), for instance. The parameters p 1 , p 2 are usually part of the trainable variables, but, as mentioned in the introduction, we shall assume that they are given (but otherwise arbitrary, e.g., p 1 and p 2 may be sampled as a random matrix and random vector respectively) for technical reasons.

In supervised learning, one seeks to tune the parameters [w, b] so that P x i (T ) approaches y i for i ∈ [N ] with respect to some loss function. To this end, the Tikhonov-regularized, empirical risk minimization problem

inf [w,b] x i (•) solves (2.4) 1 N N i=1 loss P x i (T ), y i + λ [w, b] 2 H k (0,T ;R du ) (2.8)
is considered, where λ > 0 is fixed, and x i solves (2.4) with f as in (2.5) or (2.6) (although, more general cases such as (2.7) can also be considered). Here

loss(•, •) : R m × Y → R +
is a given continuous function, which in our work we will choose on a case-by-case basis. Common examples include loss(x, y) = xy p p with p = 1, 2 for regression tasks (Y ⊂ R m ), and the cross-entropy loss

loss(x, y) = -log e xy m j=1 x j x ∈ R m , y ∈ [m]
for classification tasks. Note that (2.8) is the empirical and regularized approximation of the expected risk minimization problem:

inf [w,b] x i (•) solves (2.4) E loss (P x • (T ), •) , where E[f (•, •)] := R d ×R m f (x, y
) dρ(x, y), with x x denoting the solution to (2.4) with initial datum x. Here ρ : R d × R m → [0, 1] is an unknown joint probability distribution, from which one samples the training dataset { x i , y i } i∈ [N ] . We shall solely focus on the empirical problem in this work. By virtue of the direct method in the calculus of variations, the existence 4 of minimizers for the learning problems we consider herein can readily be shown when the appropriate parameter regularization is used (see Lemma 7.2). We put emphasis on the following remark (see also Remark 6).

Remark 1 (Sobolev regularization). We stress the consideration of a Sobolev, H 1regularization in the case of (3.2) as this is sufficient to guarantee the existence of a minimizer. Indeed, an issue could arise if solely an L 2 -regularization is used due to the specific nonlinear form of the neural ODE (3.2), which could (but might not) be an impediment for passing to the limit in the equation using only weak convergences (recall {sin(nx)} ∞ n=1 and σ(x) = max{x, 0}). This issue is specific to the continuous-time setting, as in the discrete-time thus finite dimensional optimization setting, weak and strong convergences coincide. To our knowledge, a proof of existence of a minimizer in the context of L 2 -regularized problems under System (3.2) is not known.

Empirical risk minimization

Throughout the paper we will focus on neural ODEs given by (2.4) with f as in (2.5) or (2.6). We comment on further extensions on a case-by-case basis. It will be rather convenient to work with the full stacked state trajectory

x(t) =    x 1 (t) . . . x N (t)    ,
and we recall that we make use of the convention (1.7), and we recall that d u := d×(d+1) and d x := dN . Moreover, given w ∈ R d×d and b ∈ R d , we shall write

w :=    w . . . w    ∈ R dx×dx , b :=    b . . . b    ∈ R dx . (3.1)
In view of the above discussion and noting (3.1), we will consider stacked neural ODEs in R dx such as ẋ(t) = σ(w(t)x(t) + b(t)) for t ∈ (0, T )

x(0) = x 0 ∈ R dx , (3.2) and ẋ(t) = w(t)σ(x(t)) + b(t) for t ∈ (0, T ) x(0) = x 0 ∈ R dx . (3.3)
In this section, we consider the problem of regularized empirical risk minimization. For simplicity of notation, we henceforth denote the empirical risk by

E(x) := 1 N N i=1 loss P x i , y i , (3.4)
4 Uniqueness ought not to be expected (as in most deep learning tasks) due to the inherent lack of convexity. Notwithstanding, there exist cases (e.g. two-layer perceptrons) where uniqueness can be ensured by lifting the training problem as a minimization problem over measures, which may render the problem convex (see [Chizat and Bach, 2018]).

for x ∈ R dx , where P : R d → R m and loss(•, •) ∈ C 0 (R m × Y; R + ) are given -both will change with respect to the task in question (regression, classification), as discussed in Section 2.3. For fixed λ > 0, in this section we will study the behavior when T 1 of global minimizers [w T , b T ] to the functional

J λ,T (w, b) = E(x(T )) + λ [w, b] 2 H k (0,T ;R du ) (3.5) where x ∈ C 0 ([0, T ]; R dx ) is the unique solution to either (3.3) (k = 0) or (3.2) (k = 1) corresponding to the parameters [w, b] ∈ H k (0, T ; R du ), noting (3.1).
3.1. Non-negative losses. We begin by considering the case wherein the function (x, y) → loss(P x, y) defining the empirical risk in (3.4) may attain its minimum 0. Namely, unless and until stated otherwise, we shall suppose Assumption 3. We suppose that P : R d → R m and loss ∈ C 0 (R m × Y; R + ) are such that E defined in (3.4) may attain its minimum 0. This is the case for several losses used in practice, including

• p -losses: loss(x, y) = xy p p , p ∈ {1, 2}, with P : R d → R m an affine map, or, more generally, losses which are radial functions with respect to xy (such as those induced by a distance inferred by a norm). Such modeling assumptions are typically made in the context of regression tasks ( y i ∈ Y ⊂ R m ), wherein by minimizing the empirical risk one looks to interpolate the training data by means of the projected neural ODE flow.

• The multi-class hinge loss

loss(x, y) = m j=1 j =y max {0, 1 -x y + x j } x ∈ R m , y ∈ [m],
with P : R d → R m being a matrix, which may be used for classification tasks. With regard to Assumption 3, we define the following notion of interpolation of the dataset

{ x i , y i } i∈[N ] .
Definition 3.1 (Interpolation). Let P : R d → R m be any non-zero map and let loss ∈ C 0 (R m × Y; R + ) be such that Assumption 3 is satisfied. We say that

(3.3) (resp. (3.2)) interpolates the dataset { x i , y i } i∈[N ] in some time T > 0 if there exists a time T > 0 and some parameters [w, b] ∈ L 2 (0, T ; R du ) (resp. in H 1 (0, T ; R du )) such that the unique solution x(•) to (3.3) (resp. (3.2)), noting (3.1), satisfies E(x(T )) = 0.
Note that this is a slight abuse of terminology, as precisely interpolation would rather refer to having P x i (T ) = y i for all i ∈ [N ], which is equivalent to the stated definition in the context of losses induced by a distance, but not in the context of the hinge loss. We shall, however, make use of the terminology for cohesion.

Let us also note that by means of an elementary time-scaling, if (3.3) or (3.2) interpolates the dataset in some time T > 0, it interpolates it in any time, in particular, in time 1. We will make use of this observation to simplify the subsequent presentation and analysis by assuming interpolation in time 1 without loss of generality.

We may state our main result in this context.

Theorem 3.1. Fix λ > 0. Suppose that P : R d → R m is any non-zero affine map, and suppose that loss

∈ C 0 (R m × Y; R + ) is such that Assumption 3 is satisfied. Assume that (3.3) (resp. (3.2) with σ 1-homogeneous) interpolates the dataset { x i , y i } i∈[N ] in time 1 in the sense of Definition 3.1. For any T 1 let [w T , b T ] ∈ L 2 (0, T ; R du ) (resp.
in H 1 (0, T ; R du )) be any pair of global minimizers to J λ,T defined in (3.5), and let x T (•) be the unique associated solution to (3.3) (resp. (3.2)), noting (3.1). The following properties then hold.

(i) There exists a constant

C = C { x i , y i } i∈[N ] , λ > 0 independent of T such that E(x T (T )) C T .
(ii) There exists a sequence {T n } ∞ n=1 , with T n > 0 and T n ----→ n-→∞ ∞, and some

x • ∈ R dx with E(x • ) = 0 such that, along a subsequence, x Tn (T n ) ----→ n-→∞ x • . (3.6) (iii) For any n 1, set w n (t) := T n w Tn (t T n ) for t ∈ [0, 1], b n (t) := T n b Tn (t T n ) for t ∈ [0, 1].
Then along a subsequence,

[w n , b n ] -[w * , b * ] H k (0,1;R du ) ----→ n-→∞ 0,
where

[w * , b * ] ∈ H k (0, 1; R du ) is a solution to the minimization problem inf [w,b]∈H k (0,1;R du ) x(•) solves (3.2) (resp. (3.3)) in [0,1] and E(x(1)) = 0 [w, b] 2 H k (0,1;R du ) .
Idea of proof. The proof of Theorem 3.1 may be found in Section 7.1. Let us motivate the main underlying idea.

Under the above assumptions, both (3.2) and (3.3) will be 1-homogeneous with respect to the parameters [w(t), b(t)]. Namely, both (3.2) and (3.3) (noting (3.1)) can be written as

ẋ(t) = f ([w(t), b(t)], x(t)) in (0, T ) x(0) = x 0 , (3.7) where f([αw, αb], x) = αf([w, b], x) for α > 0.
Whilst in the case of (3.3) this homogeneity property holds for any activation function σ, we require σ to be 1-homogeneous for neural networks such as (3.2). This includes rectifiers, but excludes sigmoids.

A simple computation (see Lemma 7.1) then leads to noting that, given some parameters u 1 := [w 1 , b 1 ] and the solution x 1 to

ẋ1 (t) = f w 1 (t), b 1 (t) , x 1 (t)
in (0, 1)

x 1 (0) = x 0 , (3.8) then u T (t) := 1 T u 1 ( t T ) for t ∈ [0, T ] is such that x T (t) := x 1 ( t T ) solves (3.7).
Under the interpolation assumption, we may find u 1 ∈ H k (0, T 0 ; R du ) such that the corresponding solution x 1 satisfies E(x 1 (1)) = 0, and then use the above scaling and the optimality of u T to deduce

J λ,T (u T ) E x 1 (1) + λ T u 1 2 H k (0,1;R du ) = λ T u 1 2 H k (0,1;R du )
for T 1. This will imply the decay estimate of E, and combined with some more technical compactness arguments, will yield the remaining convergence results as well.

On another hand, considering the case of (3.3) and thus k = 0, we see that

inf u T =[w T ,b T ]∈L 2 (0,T ;R du ) x T (•) solves (3.7) E(x T (T )) + λ T 0 u T (t) 2 dt = inf u T =[w T ,b T ]∈L 2 (0,T ;R du ) x T (•) solves (3.7) E(x T (T )) + λ T 1 0 T u T (sT ) 2 ds = inf u 1 =[w 1 ,b 1 ]∈L 2 (0,1;R du ) x 1 (•) solves (3.8) E(x 1 (1)) + λ T 1 0 u 1 (s) 2 ds.
This computation indicates that one may consider the behavior when T -→ ∞ for fixed λ > 0 and that when λ 0 for fixed T > 0 in the same fashion. Although this scaling is specific to the L 2 -regularization setting, it motivates completing Theorem 3.1 with the following result.

Theorem 3.2. Under the assumptions of Theorem 3.1, fix T > 0, and for any λ > 0, let [w λ , b λ ] ∈ L 2 (0, T ; R du ) (resp. H 1 (0, T ; R du )) be any pair of global minimizers to J λ,T defined in (3.5), and let x λ be the unique associated solution to (3.3) (resp. (3.2)), noting (3.1). The following properties then hold.

(i) There exists a constant C = C { x i , y i } i∈[N ] , T > 0 independent of λ > 0 such that E(x λ (T )) C λ.
(ii) There exists a sequence {λ n } ∞ n=1 , with λ n > 0 and λ n ---→ n→∞ 0, and some

x • ∈ R dx with E(x • ) = 0 such that, along a subsequence x λn (T ) ----→ n-→∞ x • .
(iii) Moreover, along a subsequence,

[w λn , b λn ] -[w * , b * ] H k (0,T ;R du ) ----→ n-→∞ 0, where [w * , b * ] ∈ H k (0, T ; R du ) is a solution to the minimization problem inf [w,b]∈H k (0,T ;R du ) x(•) solves (3.2) (resp. (3.3)) and E(x(T )) = 0 [w, b] 2 H k (0,T ;R du ) .
Remark 2 (Homogeneity).

• The impediment that appears in the results presented above when considering neural ODEs of the form

ẋ(t) = w 1 (t)σ w 2 (t)x(t) + b 2 (t) + b 1 (t) in (0, T ) x(0) = x 0 , (3.9)
where w 1 (t) ∈ R d×d hid and w 2 ∈ R d hid ×d , with w 1 (t) and w 2 (t) defined after stacking the states {x i (t)} i∈[N ] within x(t) following convention (1.7), is the lack of homogeneity (and thus scaling) of the dynamics with respect to the parameters. Let us elaborate by focusing on b 2 ≡ b 1 ≡ 0 for simplicity. Let x T denote the solution to (3.9) with parameters

[w 1 T , w 2 T ]. We set x 1 (t) := x T (tT ) for t ∈ [0, 1],
and we see that, assuming σ is α-homogeneous with α > 0, x 1 solves (3.9) on [0, 1] with parameters

w 1 1 , w 2 1 := T p w 1 T (tT ), T q /α w 2 T (tT )
for t ∈ (0, 1), with p + q = 1, noting (3.1). Solely looking at the squared L 2 -norms of the parameters, we see that

T 0 w 1 T (t) 2 dt + T 0 w 2 T (t) 2 dt = T 1 0 w 1 T (sT ) 2 ds + T 1 0 w 2 T (sT ) 2 ds = T 1-2p 1 0 w 1 1 (s) 2 ds + T 1-2q /α 1 0 w 2 1 (s) 2 ds.
We see that, in order to have a scaling factor which decays when T -→ ∞, a property which is a cornerstone of our proof, we would simultaneously need 1 -2p < 0 and 1 -2q α < 0, i.e. p > 1 2 and q > α 2 . Due to the constraint p + q = 1, this would entail α < 1, which means that σ cannot be globally Lipschitz continuous.

• Let us, in line with what is discussed just above, check what is the needed degree of homogeneity of σ in the case of System (3.2) to ensure that the results of Theorem 3.1 hold. Let us thus assume that σ is α-homogeneous with α > 0. Then, setting x 1 (t

) := x T (tT ) for t ∈ [0, 1] where x T solves (3.2) with parameters [w T , b T ], we see that x 1 (t) solves (3.2) on [0, 1] with parameters [w 1 (t), b 1 (t)] = T 1 /α w T (tT ), T 1 /α b T (tT ) for t ∈ (0, 1).
Then, looking at the L 2 -norms of the parameters, we see that

T 0 w T (t) 2 dt + T 0 b T (t) 2 dt = T 1 0 w T (sT ) 2 ds + T T 0 b T (sT ) 2 ds = T 1-2 /α 1 0 w 1 (s) 2 ds + T 1-2 /α 1 0 b 1 (s) 2 ds,
and thus, our methodology really needs α ∈ [1, 2). Note that this brief computation also indicates the precise proximity of E(x T (T )) to 0 for σ is αhomogeneous, which is of the order of O T 1-2 /α .

Remark 3 (The output layer P ). We note that the output layer parameters given by the affine map P are fixed, but in general arbitrary. This is due to the fact that if we were to optimize P as well, we would have to ensure that the optimal P is bounded with respect to the limiting hyper-parameter (λ or T ). This in turn could perhaps be ensured if we were to regularize the output layer as well, but would, in turn, be an impediment to the scaling arguments we use in all proofs since then the parameter regularization norm would not scale polynomially with T .

3.2. Positive losses. We henceforth consider the standard setting of classification tasks, wherein the labels y i take values in a set of m 2 classes -unless stated otherwise, we consider y i ∈ [m] for all i ∈ [N ]. We will focus on the cross-entropy loss in (3.4), which we recall, when evaluated along the neural ODE output, reads

loss(P x i (T ), y i ) := -log e P x i (T ) y i m j=1 e P x i (T ) j .
(3.10)

Here P : R d -→ R m is an affine map, with more assumptions made precise later on.

An important feature of the cross-entropy loss is the fact that it is not coercive with respect to the first variable -namely, as P x i (T ) y i goes to infinity, the loss goes to zero. In particular, the infimum 0 is never attained. This is in line with intuition regarding the classification task, as one looks to separate the features with respect to their individual class in the label space R m .

The problem consisting of classifying a given dataset is closely tied to the following rather intuitive notion of separability, which we will need in subsequent results. Definition 3.2 (Separability). Let P : R d -→ R m be any non-zero affine map. We say that

(3.3) (resp. (3.2)) separates the dataset { x i , y i } i∈[N ] with respect to P if there exists a time T > 0 and some parameters [w, b] ∈ L 2 (0, T ; R du ) (resp. in H 1 (0, T ; R du )) such that the unique solution x(•) to (3.3) (resp. (3.2)) satisfies P x i (T ) y i -max j∈[N ] j = y i P x i (T ) j > 0 for all i ∈ [N ].
In other words, a parametrized neural ODE flow separates the given dataset if the corresponding margin γ [w,b] , defined as

γ [w,b] := min i∈[N ]   P x i (T ) y i -max j∈[N ] j = y i P x i (T ) j    (3.11)
is positive. We may now state our main result in the classification context, which entails a quantitative rate of decay as T -→ ∞ of the training error with cross-entropy loss for ReLU activated neural ODEs.

Theorem 3.3. Let { x i , y i } i∈[N ] be a given dataset with x i ∈ R d and y i ∈ [m]. Let λ > 0 be fixed, and let Q : R dx -→ R d be a non-zero affine map such that Q x i 0 for i ∈ [N ].
Set (recall convention (1.7))

x 0 i := Q x i for i ∈ [N ],
and let P ∈ R m×d be any non-zero matrix such that System (3.2), with σ(x) = max{x, 0}, separates the dataset x 0 i , y i i∈[N ] with respect to P in some time T 0 > 0 as per Definition 3.2, and let γ denote the associated margin as defined in (3.11). For any T T 0 , let [w T , b T ] ∈ H 1 (0, T ; R du ) be any pair of global minimizers to J λ,T defined in (3.5)-(3.10), and let x T (•) be the associated unique solution to (3.2) with σ(x) = max{x, 0}.

Then, there exists a constant

C = C { x i , y i } i∈[N ] , λ > 0 independent of T > 0 such that E(x T (T )) log 1 + (m -1)e -γ e T α 2 + C T 2α-1
(3.12)

holds for any α ∈ 0, 1 2 . By using a Taylor expansion, one sees that the first term in the upper bound is negligible when T 1. Thus, the training error E(x T (T )) is at most of the order O T 2α-1 .

We note that the above theorem is very specific to neural ODEs of the form (3.2) with ReLU activations, and the specific form of the cross-entropy loss, from which the first term in the estimate (3.12) is derived. This is due to the proof strategy, which relies on using the positivity of the right hand side to, in some sense, obtain a linear equation for the projected output features for some auxiliary parameters constructed within the proof, and thus have an explicit solution for these parameters of the form ∼ e t . This stimulates the appearance of the second exponential within the log in (3.12).

Much like what we observed in the setting of losses which attain their minimum, we can expect to link the limit as T goes to infinity with the convergence of the regularization path, namely the limit as λ 0. This is depicted in the following theorem.

Theorem 3.4. Under the assumptions of Theorem 3.3, fix T T 0 and for any λ > 0 let [w λ , b λ ] ∈ H 1 (0, T ; R du ) be any pair of global minimizers to J λ,T defined in (3.5)-(3.10), and let x λ (•) be the associated unique solution to (3.2) with σ(x) = max{x, 0}. Then, there exists a constant

C = C { x i , y i } i∈[N ] , T > 0 independent of λ > 0 such that E(x λ (T )) log 1 + (m -1)e -γ e λ -α 2 + C λ -2α+1
holds for any α ∈ 0, 1 2 . Remark 4 (Absence of margin convergence). We note that, unlike Theorem 3.1 (and Theorem 3.2), in Theorem 3.3 (and Theorem 3.4) we do not provide any result regarding the behavior of the optimal parameters [w

T , b T ] as T -→ ∞ (resp. [w λ , b λ ] as λ 0)
. Let us elaborate on this absence in the case λ 0.

• In [START_REF] Wei | Regularization matters: Generalization and optimization of neural nets vs their induced kernel[END_REF] (see also [START_REF] Rosset | Margin maximizing loss functions[END_REF]), the authors show that for the perceptron without bias 5 Φ(x, u) := w 2 σ(w 1 x), where σ is, say, 1-homogeneous, u = w 1 , w 2 with w 2 ∈ R m×d hid and w 1 ∈ R d hid ×d , any global minimizer u λ to

J λ (u) := 1 N N i=1 -log e Φ( x i ,u) y i m j=1 e Φ( x i ,u) j + λ u 2 is such that γ u λ ---→ λ 0 γ * ,
where u λ = uλ / uλ , γ u λ is the margin of Φ(•, u λ ) as defined in (3.11), and γ * > 0 is the max-margin:

γ * := max u 1 γ u .
5 In fact, the result holds for a multi-layer perceptron of arbitrary but fixed depth -we concentrate on this simpler scenario for presentational purposes.

A set of parameters u * maximizing the margin are then called a max-margin solution (here we focus on 2 -margins; the 1 -case has also been studied in the machine literature). The proof of this fact relies crucially on the specific form of the perceptron, as the map u -→ Φ(x, u) is positively homogeneous (of degree 2) for any x ∈ R d . By using the exponential scaling of the cross entropy, J λ (u λ ) can be lower bounded roughly by e -u λ 2 γ λ , and also has an upper bound that scales with e -u λ 2 γ * . Again making use of the fact that u -→ Φ(x, u) is positively homogeneous, one can ensure that u λ -→ ∞ as λ 0, and then take u λ large enough so that γ *γ λ vanishes.

In the context of neural ODEs such as (3.2) (and the corresponding ResNet, for that matter), to replicate such ideas, one would need to ensure that the map [w, b] -→ P x(T ), where x solves (3.2), is positively homogeneous (of some degree). In fact, let us see under which conditions the following relaxed problem, can be solved. We seek a non-decreasing continuous function f : (0, ∞) → (0, ∞) such that

x α = f (α)x
for all α > 0,

where
(t) = f (α) ẋ(t) = f (α)σ(w(t)x(t)).
Using the fact that ẋα (t) solves (3.2) with parameters [αw, 0], we moreover see that σ(αw(t)f (α)x(t)) = f (α)σ(w(t)x(t)).

Now we see that the above relation translates to having

σ(αf (α)s) = f (α)σ(s) for all α > 0, s ∈ R. (3.13)
This implies σ(0) = 0, and would also be an obstruction to the Lipschitz continuity of σ at x = 0 (indeed, fix s ∈ R such that σ(s) = 0, and let α 0). Thus, the only increasing function f : (0, ∞) → (0, ∞) such that (3.13) holds and σ ∈ Lip(R) is the zero function.

• The issue presented just above does not occur when the only trainable parameter is the additive bias b(t). Indeed, if for instance

w(t) = w • = diag(w • ) is fixed, x(t) solves (3.2) with parameters [w • , b(t)] and x α (t) solves (3.2) with parameters [w • , αb(t)]
, then one readily sees that P x α,i (t) = αP x i (t) for all α > 0 whenever P x i = 0 for i ∈ [N ] (this is the case when augmenting the data x i and projecting orthogonally onto the added euclidean space via P ). Hence, one may show, by adapting the techniques of [START_REF] Wei | Regularization matters: Generalization and optimization of neural nets vs their induced kernel[END_REF], that the margin γ [w•,b T ] corresponding to normalized global minimizers b T = bT / bT , converges to the max-margin

γ * := sup b H 1 (0,1;R du ) 1 γ [w•,b]
as T -→ ∞ (with an analog result when λ 0 and T is fixed). We omit the statement and proof of this result due to the possible limitations of not optimizing the weights in practical contexts.

In view of this discussion, one sees that in order to deduce what the limit for the margin in the fully trained neural ODE (and ResNet) setting would be, one would likely have to conceive a different strategy to that of [START_REF] Rosset | Margin maximizing loss functions[END_REF][START_REF] Wei | Regularization matters: Generalization and optimization of neural nets vs their induced kernel[END_REF] which relies on the homogeneity of the parameter to output map.

Remark 5 (Different losses). We note that in the context of binary classification, one could consider y i ∈ {-1, 1} and train with the logistic loss in (3.4) loss(P x i (T ), y i ) := log 1 + e -y i P x i (T ) , where P : R d → R is an affine map. The statements of Theorem 3.3 and Theorem 3.4 also hold in the context of the logistic loss defined above by a straightforward repetition of our proofs.

Another class of losses when y i ∈ [m] is the purely exponential loss loss(P x i (T ), y i ) := e -y i P x i (T ) .

Note that, again, by a slight adaptation of our techniques, one can show an analog variant of Theorem 3.3 (and Theorem 3.4) with an estimate of the form

E(x T (T )) e -γe T α 2 + C T 2α-1 .

Augmented empirical risk minimization

We are now interested in seeing whether one can obtain better quantitative estimates for the decay of the training error E to 0 with respect to the time horizon (∼ number of layers) T > 0 -namely, improve the O 1 T -rate of convergence of the training error to 0 manifested in Theorem 3.1 and Theorem 3.3. We provide a proof of exponential stability in the setting of an augmented supervised learning problem, concentrating solely on L 2parameter regularization (and hence we shall assume the existence of a minimizer in the setting of (3.2)). We shall also consider a BV-parameter regularization, which suffices to ensure the existence of minimizers, but only provide a polynomial convergence of the averages of the training error and optimal parameters to zero as T -→ ∞, whilst ensuring a uniform bound on the total variation of the optimal parameters. Unless stated otherwise, we will henceforth solely concentrate on empirical risks

E(x) := 1 N N i=1 loss(P x i , y i ) for x ∈ R dx , (4.1)
induced by a loss satisfying the following condition.

Assumption 4. We assume that Y ⊂ R m and loss ∈ C 0 (R m × Y; R + ) is such that there exist constants c > 0 (possibly depending on m) and α > 0 such that

loss(x, y) c x -y α 2 for x, y ∈ R m .
Typical examples of such losses are loss(x, y) = xy p p for p 1 and Y ⊂ R m ; note however that this assumption excludes the log-loss and the cross-entropy loss. In (4.1), P ∈ Lip(R d ; R m ) is any given surjective and non-zero map, which, in the context of regression tasks, is additionally an affine map, while in the context of binary classification, e.g. y i ∈ Y = {-1, 1}, may be an affine map composed with a thresholding nonlinearity with range [-1, 1].

Inspired from insights in optimal control theory, for fixed λ > 0, we will study the behavior when T 1 of global minimizers to the functional

J T (w, b) := E(x(T )) + 1 N T 0 x(t) -x 2 dt + λ [w, b] 2 L 2 (0,T ;R du ) , (4.2)
with E as in (4.1)-Assumption 4, and where x i ∈ P -1 ({ y i }) for all i ∈ [N ] are arbitrary but fixed. We note that, contrary to the case where we minimize the training error at the final time T , here the same scaling does not appear, which allows us to deduce an equivalence with λ -→ 0. Hence, we will solely be interested in the behavior when T 1. In fact, we ought to expect a result of slightly different nature to Theorem 3.1. This is because the trajectory tracking term introduces a stronger time-scale in the behavior of the optimization problem as T 1. To see this, consider the neural ODE (3.3) for simplicity, and note that

inf u T ∈L 2 (0,T ;R du ) x T solves (3.7) E(x T (T )) + 1 N T 0 x T (t) -x 2 dt + λ T 0 u T (t) 2 dt = inf u T ∈L 2 (0,T ;R du ) x T solves (3.7) E(x T (T )) + T N 1 0 x T s T -x 2 ds + λ T 1 0 T u T (sT ) 2 ds = inf u 1 ∈L 2 (0,1;R du ) x 1 solves (3.8) E x 1 (1) + T N 1 0 x 1 (s) -x 2 ds + λ T 1 0 u 1 (s) 2 ds. (4.3)
We see that, unlike Theorem 3.1, the trajectory tracking term in (4.3) carries significance when T 1.

We will require the following controllability notion, which is rather natural in the context of the result that follows.

Definition 4.1 (Controllability with linear cost). We say that • System (3.3) is controllable with linear cost if for any T > 0 and r > 0 there exists a constant C(T, r) > 0 such that for any x 0 ∈ R dx and x 1 ∈ R dx satisfying x 0x 1 r, there exists a pair of parameters [w, b] ∈ L 2 (0, T ; R du ) for which the corresponding unique solution x(•) to (3.3), noting (3.1), satisfies x(T ) = x 1 , and

[w, b] L 2 (0,T ;R du ) C(T, r) x 0 -x 1 (4.4)
holds.

• System (3.2) is controllable with linear cost if for any T > 0 and r > 0 there exists a constant C(T, r) > 0 such that for any x 0 ∈ R dx and x 1 ∈ R dx satisfying x 0x 1 r, there exists a pair of parameters [w, b] ∈ C 0 ([0, T ]; R du ) ∩ BV(0, T ; R du ) for which the corresponding unique solution x(•) to (3.2), noting (3.1), satisfies x(T ) = x 1 , and

[w, b] L ∞ (0,T ;R du ) + [w, b] BV(0,T ;R du ) C(T, r) x 0 -x 1 (4.5)
holds.

We refer to Theorem 5.1 for further analysis and comments regarding Definition 4.1, in particular regarding the amplitude estimates (4.4) -(4.5). Note that (4.5) is a weaker condition than assuming an H 1 (0, T ; R du )-estimate due to the Sobolev embedding.

On another hand, due to the homogeneity of the dynamics, one readily sees (as per Lemma 7.1) that the controllability assumption holds in any time T > 0 if it holds in some time T 0 > 0.

We refer to (4.4) -(4.5) as linear estimates on the cost of [w, b] as, in addition to the linear growth with respect to x 0x 1 such estimates are a hallmark of linear controlled dynamical systems of the form ẋ(t) = Ax(t) + Bu(t), where u(t) ∈ R m is the control, x(t) ∈ R n is the state, with A ∈ R n×n and B ∈ R n×m being given matrices (see [START_REF] Zuazua | Controllability and observability of partial differential equations: some results and open problems[END_REF]).

We are in a position to state our main result in the context of the augmented supervised learning problem consisting of minimizing (4.2).

Theorem 4.1 (Exponential stability). Fix λ > 0, let P ∈ Lip(R d ; R m ) be any given non-zero surjective map, and let x ∈ R dx with x i ∈ P -1 ({ y i }) for i ∈ [N ] be arbitrary but fixed. Suppose that system (3.3) (resp. (3.2) with σ 1-homogeneous) is controllable with linear cost in the sense of Definition 4.1. Then, there exists T * > 0 and constants

C = C { x i , y i } i∈[N ] , λ, N, α, P > 0 and µ = µ { x i , y i } i∈[N ] , λ, N, α > 0 such that for any T T * , any pair of parameters [w T , b T ] ∈ L 2 (0, T ; R du ) minimizing J T defined in (4.2), and the corresponding unique solution x T (•) to (3.3) (resp. (3.2)) satisfy w T (t) + b T (t) C e -µt for a.e. t ∈ [0, T ] and E(x T (t)) + x T (t) -x C e -µt for all t ∈ [0, T ].
The convergence/stability rate entailed by Theorem 4.1 is not only noticeably stronger than those in Theorem 3.1 and Theorem 3.3, but the estimate holds in any time t ∈ [0, T ] (i.e., at every layer when viewed from the ResNet perspective) and not only for the output features.

Theorem 4.1 can be proven by following the framework presented 6 in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF]. We provide most of the arguments in Section 7.4 due to technical changes in the underlying model.

Remark 6 (Existence of minimizers, BV-regularization). In Theorem 4.1 we have subjacently assumed the existence of a minimizer of (4.2) for the neural ODE (3.2). By solely regularizing the L 2 -norm of the parameters, we are, a priori, not aware how to ensure the strong convergence in L 1 of a minimizing sequence of parameters, which would suffice for applying the direct method in the calculus of variations. On the other hand, using a Sobolev H 1 -regularization would entail that optimizable parameters are continuous (by the Sobolev embedding) and thus render our specific proof strategy incompatible, as we rely on constructing discontinuous suboptimal parameters by zero extensions 6 As discussed in the introduction, Theorem 4.1 is a specific manifestation of the so-called turnpike property, a paradigm dating back to the works of John von Neumann [START_REF] Von Neumann | A model of general economic equilibrium[END_REF]), and works in economics by Paul Samuelson et al. [START_REF] Dorfman | Linear Programming and Economic Analysis[END_REF]). A turnpike theory under smallness assumptions on the inputs and/or targets, combining the Pontryagin Maximum Principle, linearization arguments and precise estimates on Riccati equations, and covering a wide variety of nonlinear optimal control problems is developed in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] -with extensions to Lipschitz nonlinearities and avoiding smallness conditions found in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF].

In view of this, in the setting of (3.2) we may also make use of a BV-regularization rather than H 1 , as while either of these regularizations suffice to guarantee compactness of the flow, which is sufficient to ensure the existence of a minimizer, the BVregularization allows for discontinuous parameters: one would then rather consider

J T (w, b) := E(x(T )) + 1 N T 0 x(t) -x 2 dt + λ [w, b] 2 BV(0,T ;R du ) . (4.6)
Taking the above remark into consideration, we are able to prove the following partial stability and convergence result in the setting of BV-regularized parameters. We focus on (3.2), since (3.3) can be covered by means of L 2 -regularization.

Theorem 4.2 (Convergence of averages, stability). Fix λ > 0, let P ∈ Lip(R d ; R m ) be any given non-zero surjective map, and let x ∈ R dx with

x i ∈ P -1 ({ y i }) for i ∈ [N ]
be arbitrary but fixed. Suppose that (3.2), with σ 1-homogeneous, is controllable with linear cost in the sense of Definition 4.1. The following facts then hold.

(i) (Stability). There exists a constant

C = C { x i , y i } i∈[N ] ,
x 0x , λ, N, P > 0 such that for any T 1, any pair of parameters [w T , b T ] ∈ BV(0, T ; R du ) minimizing J T defined in (4.6), and the corresponding unique solution x T (•) to

(3.2) satisfy E(x T (t)) + x T (t) -x + D [w T , b T ] (0, T ) C for all t ∈ [0, T ], where |D [w T , b T ]| (0, T ) denotes the total variation of [w T , b T ].
(ii) (Convergence of averages). Any pair of parameters [w T , b T ] ∈ BV(0, T ; R du ) minimizing J T defined in (4.6), and the corresponding unique solution x T (•) to (3.2) also satisfy

1 T T 0 E(x T (t)) dt + 1 T T 0 x(t) -x dt + 1 T T 0 [w T (t), b T (t) dt ----→ T -→∞ 0.
The convergence of averages are very much related to the so-called integral turnpike property, which has been studied in some optimal control contexts (see e.g. [START_REF] Trélat | Integral and measure-turnpike properties for infinitedimensional optimal control systems[END_REF]). Note that, in addition to these convergences, the stability estimates of item (i) entail that the oscillations of the training error may be controlled, in any time t, uniformly with respect to the time horizon T . Moreover, the total variation of the optimal parameters is also uniformly bounded with respect to T . We are unable to provide exponential stability estimates as in Theorem 4.1 due to some constructions specific to the proof of Theorem 4.1, wherein suboptimal parameters with jumps at time instances are constructed. While the L 2 -norm does not see these jump singularities, the BV-norm does and renders our strategy incompatible.

Remark 7 (Extensions). We only stated Theorem 4.1 for neural ODEs of the form (3.3), or (3.2) with σ 1-homogeneous. This is solely to guarantee the exponential stability estimate of the optimal parameters, for which our proof requires using the underlying scaling endowed by the homogeneity of the dynamics. But in fact, the exponential stability estimate of the training error E(x T (T )) and a uniform-in-T bound of the optimal parameters can be shown for more complicated neural ODE dynamics such as (3.9), solely by a small adaptation of the proof.

Remark 8 (Dependence on N ). All of the results presented in what precedes hold for a fixed but otherwise arbitrary number of data samples N . In Theorem 4.1, one also notes an explicit dependence of the constants C > 0 and µ > 0 on N -in fact, both constants might have a tendency to depend in an exponential manner with respect to N due to the subjacent application of a Grönwall inequality for the stacked neural ODE system, for which the parameters are pasted N times. The dependence on N provided in our proof might not be sufficient or sharp for studying a possible large data limit (e.g., via a law of large numbers argument of some kind) -we leave this issue open.

Example 4.1 (A numerical experiment). In Figure 1 -Figure 2, we depict 7 a manifestation of the exponential stability insinuated by Theorem 4.1 on a toy binary classification task (namely y i ∈ {-1, 1}) with N = 2400 training samples and 600 test samples, where P (•) = tanh(p 1 • +p 2 ) with p 1 ∈ R 1×2 and p 2 ∈ R randomly sampled from a uniform distribution on [0, 10]. Note that, while in theory tanh : R → [-1, 1] is not surjective (it is only bijective onto (-1, 1)), we use it as a thresholding nonlinearity for simplicity, as of course numerically the lack of surjectivity does not appear due to floating point accuracy.

To discretize the full continuous-time optimization problem, we use direct shooting, which is a first discretize then optimize approach. We consider the neural ODE (3.3) with σ(x) = tanh(x) (we use the ResNet (2.3)), with T = 15 (and thus 15 layers) and λ = 0.01. Finally, we discretize the integrals using an elementary trapezoidal quadrature. We note that the learned flow has a distinctly simple variation in Figure 2, and, albeit on a toy task, we observe satisfactory generalization properties in Figure 3.

E(x T (t)) (left) over t ∈ [0, T ].
We observe that, after a finite time, the training error and trajectory remain at a steady configuration, so further times/layers could be discarded from training.

7 Software experiments were done using PyTorch [START_REF] Paszke | Automatic differentiation in PyTorch[END_REF] (and may be found at https://github.com/borjanG/dynamical.systems), using the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with learning rate equal to 10 -3 . Experiments were conducted on a personal MacBook Pro laptop (2.4 GHz Quad-Core Intel Core i5, 16GB RAM, Intel Iris Plus Graphics 1536 MB). we have restricted our study to a trajectory tracking term consisting of the squared L 2 (0, T ; R dx )-norm, even-though the final cost E(x T (T )) allows us to address both classification and regression tasks. However, having to look for targets x in the preimage of the labels y i by P for any general task may not scale computationally and might bias the prediction.

Interestingly enough, at least numerically (any analytical result remains an open problem), we observe that the stabilization phenomenon persists (although, in theory, perhaps not with the same rate) when the term x(t)-x 2 is replaced by the training error E(x(t)) with a general and possibly non-coercive loss, for instance, the cross-entropy loss on a multi-label classification tasks as seen in Figure 7 & Figure 13. In fact, one could stipulate that this stabilization phenomenon (be it exponential or not, but in any t ∈ [0, T ] rather than just for the output features at time T ) holds for global minimizers of functionals of the form

J T (w, b) := T 0 E(x(t)) dt + λ [w, b] 2 H(0,T ;R du ) , (4.7)
with E as in (3.4) and loss being continuous and nonnegative, but otherwise arbitrary; H is for instance L 2 or BV. We perform several numerical experiments to justify 8 this claim (see Example 4.3 -Example 4.6). Note that in these experiments, even-though the loss is taken as cross-entropy and is thus not coercive, in addition to a stability property for E(x(t)) to 0 we also see (e.g. in Figure 7) that the trajectories x(t) and features {P x i (t)} i∈[N ] stabilize towards some targets in sufficiently large time. However, due to the fact that E does not attain its minimizer, it is a priori not clear how one may characterize the targets to which x(t) and features {P x i (t)} i∈[N ] stabilize. We refer to [START_REF] Yagüe | Sparse approximation in learning via neural ODEs[END_REF] for a result illustrating a similar stability phenomenon in the setting of (4.7) with L 1 -parameter regularization.

Example 4.2 (Concentric spheres). For comparison purposes, we consider an identical dataset setting to that of Figure 1 -Figure 3. We consider the neural ODE (3.3) with σ ≡ tanh, cross-entropy loss and λ = 0.01, with the output layer having the form P x = p 1 x + p 2 , with p 1 ∈ R 2×2 , p 2 ∈ R 2 both being part of the trainable parameters. We visualize the output of the experiments in Figure 4 -Figure 6 below. We do not insinuate that our numerical experiments are comparable with state of the art configurations, as we only look to motivate and illustrate the mathematical phenomena studied in what precedes. Indeed, we generally make use of a forward Euler scheme compared to more advanced adaptive schemes, as used for instance in [START_REF] Chen | Neural ordinary differential equations[END_REF]. Our experiments should rather be seen as proof of concept. The trained classifier on [-2, 2] 2 and its evaluation on a batch of the test data. The shape of the data is captured accurately. data is used, and the trajectories evolve in the ambient dimension d = 2. The output layer is parametrized by P x = p 1 x + p 2 , where p 1 ∈ R 3×2 , p 2 ∈ R 3 are part of the trainable parameters. We display the results of the experiments in Figure 7 -Figure 9.

Example 4.4 (XOR). We consider a canonical binary classification task ( y i ∈ [2])the XOR dataset (Figure 11) consisting of N = 3200 training samples and 800 test samples. To further illustrate the genericity of the stability phenomenon, we now consider the neural ODE (3.9), with T = 15 (we use a forward Euler scheme to obtain a corresponding ResNet with fixed time-step equal to 1), where σ ≡ tanh and d hid = 3. We again consider cross-entropy loss in the empirical risk E defined in (4.7) and λ = 0.01, and we focus solely on L 2 -parameter regularization (for simplicity). No augmentation of the initial data is used, and the trajectories evolve in the ambient dimension d = 2.

The output layer is defined as P x = p 1 x + p 2 , with p 1 ∈ R 2×2 , p 2 ∈ R 2 both being part of the trainable parameters. We display the results in Figure 10 -Figure 12.

Example 4.5 (MNIST). We now show that the stabilization phenomenon may also be observed on more realistic datasets such as MNIST [LeCun et al., 2010]. MNIST is a dataset consisting of handwritten digits ranging from 0 to 9, with a training set consisting of 60000 samples, and a test set consisting of 10000 samples. Each input sample x i is a The trained classifier on [-2.5, 2.5] 2 and its evaluation on a batch of the test data. The shape of the dataset is captured accurately. We depict the evolution of two individual samples x i (t) 2 R 784 mapped onto a 28 ⇥ 28 grid. We see that each trajectory stabilizes to some stationary configuration, and the trained model tends to compress the input digit samples ahead of classifying them via the softmax applied to P x i (t) 2 R 10 .

regime (E(x(T )) = 0 with E given in (3.4)) when T increases. It is thus of interest to also illuminate some of the properties of the parameters which allow the trajectory to reach a minimizer of the empirical risk E, and to see whether such parameters indeed exist.

By means of an elementary Grönwall argument, we can show the following illustrative result, which stipulates a lower bound for the amplitude of the weights w in terms of We depict the evolution of two individual samples x i (t) ∈ R 784 mapped onto a 28 × 28 grid. We see that each trajectory stabilizes to some stationary configuration. The trained model tends to "diffuse" (in a caloric sense) the input signal ahead of classifying via the softmax applied to P x i (t) ∈ R 10 .

The interpolation regime

The majority of our results stated in the preceding sections stipulate whether and how the output x(T ) of the neural ODE trajectory approaches the so-called interpolation regime (E(x(T )) = 0 with E given in (3.4)) when T increases. It is thus of interest to also illuminate some of the properties of the parameters which allow the trajectory to reach a minimizer of the empirical risk E, and to see whether such parameters indeed exist.

By means of an elementary Grönwall argument, we can show the following illustrative result, which stipulates a lower bound for the amplitude of the weights w in terms of the dispersion or concentration of the input data.

Proposition 5.1. Let P : R d → R m be surjective, and let T > 0. Assume that for some parameters [w, b] ∈ L 1 (0, T ; R du ) the solution x ∈ C 0 ([0, T ]; R dx ) to either (3.3) The validation error and accuracy over training epochs (experiments repeated 10 times); as anticipated, generalization is not as good as for the simpler MNIST dataset. The lower accuracy with respect to state of the art configurations could also be due to the fact that we do not make use of convolutional layers.

or (3.2) satisfies P x i (T ) = y i for all i ∈ [N ].
(5.1)

Then w L 1 (0,T ;R du ) C(σ) max (i,j)∈[N ] 2 i =j inf x 1 i ∈P -1 ({ y i }) x 1 j ∈P -1 ({ y j }) log   x 1 i -x 1 j x 0 i -x 0 j   , (5.2)
where C(σ) > 0 is the Lipschitz constant of σ ∈ Lip(R).

By virtue of Cauchy-Schwarz, (5.2) clearly implies

w L 2 (0,T ;R du ) C(σ) √ T max (i,j)∈[N ] 2 i =j inf x 1 i ∈P -1 ({ y i }) x 1 j ∈P -1 ({ y j }) log   x 1 i -x 1 j x 0 i -x 0 j   .
This observation leads us to stipulate the interest of L 2 -regularization and increasing T : should the training data contains inputs which are very concentrated in the ambient space, parameters in the interpolation regime ((5.1) is equivalent to E(x(T )) = 0 whenever E attains its minimum) will dissipate in L 2 -norm over longer time horizons. On another hand, the difference quotient

κ x 0 i , x 1 i := x 1 i -x 1 j x 0 i -x 0 j
can roughly be seen as an indicator of the variations of the flow, and thus the weight matrix would control the latter's simplicity. Indeed, denoting by Φ T : R d → R d the map defined by Φ T (x 0 i ) = x i (T ), where x i (t) solves either (3.3) or (3.2), one can roughly stipulate that κ is an approximation ∇Φ T (evaluated at some intermediate point -here ∇Φ T denotes the Jacobian matrix of the flow map Φ T ). This can be argued already by using, for instance, the Cauchy mean-value theorem. Let us provide further detail and focus on (3.3) for simplicity, and assume σ ∈ C 1 (R) ∩ Lip(R). We may linearize (3.3) with respect to the initial datum

x 0 i to obtain                żi (t) = w(t)      σ x 0 i,1 z i,1 (t) . . . σ x 0 i,d z i,d (t)      := w(t)z i (t) for t ∈ (0, T ) z i (0) = z 0 i for i ∈ [N ]. One then sees that exp T 0 w(t) dt = ∇Φ T (x 0 i ) for i ∈ [N ].
One may thus see T 0 w(t) dt as a "multi-dimensional" logarithm of ∇Φ T (x 0 i ). Thus, the Jacobian of the flow, which is an indicator of its variations and thus of its simplicity, can be measured by the weight matrix. 5.1. On Definition 4.1. To complete this section, we state the following interpolation result, which includes an estimate on the parameters with respect to the distance of the target and the initial datum, which somewhat enhances the validity of the assumption we make in Theorem 3.1. While such an estimate is standard in the setting of linear models, it is not provided by sufficient controllability conditions for nonlinear systems such as the Chow-Rashevski theorem [Coron, 2007, Chapter 3, Section 3.3].

Theorem 5.1. Let T > 0 and assume that N d. Let x 1 ∈ R dx be given, and assume that the activation function σ ∈ C 1 (R) ∩ Lip(R) is such that σ x 1 1 , . . . , σ x 1 i , . . . , σ x 1 N is a system of linearly independent vectors in R d . Then, there exist universal constants r > 0 and C > 0 such that for any datum x 0 ∈ R dx satisfying x 0x 1 r, there exists a weight matrix w ∈ L ∞ (0, T ; R d×d ) such that the unique solution x(•) to

ẋ(t) = w(t)σ(x(t)) in (0, T ) x(0) = x 0 , satisfies x(T ) = x 1 ,
and the following estimate holds

w L ∞ (0,T ;R d×d ) C T x 0 -x 1 .
Remark 9.

• For simplicity of presentation, we have not exhibited the bias parameter b(t). One can readily check that, in the presence of this additional parameter, the assumption N d can be relaxed to N d + 1.

• The case N > d + 1 may be treated by, for instance, appending additional features (e.g. zeros as in [START_REF] Dupont | Augmented Neural ODEs[END_REF]) to the input data x 0 i for i ∈ [N ] as to guarantee that the augmented datum is of dimension d aug N .

In the discrete-time context of neural networks such as (2.1) or (2.3), the property analog to Definition 3.1 is also well explored in the literature, and is commonly called finite sample expressivity [START_REF] Zhang | Understanding deep learning requires rethinking generalization[END_REF]). An additional interest is that of estimating the number of parameters -referred to as the memorization capacity -needed to manifest this property. For further results in this direction we refer the reader to [START_REF] Yun | Small ReLU networks are powerful memorizers: a tight analysis of memorization capacity[END_REF]Bubeck et al., 2020a,b] and the references therein.

In the ODE context, the property of finite sample expressivity finds its analog in the ensemble or simultaneous controllability, wherein one requires only 1 pair of parameters/controls to steer N trajectories of the same system to N prescribed targets -this is the property we show in Theorem 5.1. The literature on such controllability results of neural ODEs, mostly relying on geometrical techniques such as Lie brackets techniques (see [Coron, 2007, Chapter 3, Section 3.3]), under specific constraints on the activations function, is vast (see e.g. [START_REF] Cuchiero | Deep neural networks, generic universal interpolation, and controlled ODEs[END_REF][START_REF] Agrachev | Control on the manifolds of mappings as a setting for deep learning[END_REF][START_REF] Tabuada | Universal approximation power of deep neural networks via nonlinear control theory[END_REF]). We refer to [START_REF] Ruiz-Balet | Neural ODE control for classification, approximation and transport[END_REF] for further results in this direction.

Continuous space-time neural networks

We now come back to the scheme (2.3) defining a ResNet with N layers 2 layers. Whilst such networks are widely used in practice, in the discrete-time context, they do not take into account variations of the dimensions of the weights and states over layers. Such variations may arise when considering convolutional and/or pooling layers, which are ubiquitous in tasks in computer vision. In such tasks, it is moreover of interest to view the data itself as being continuum objects.

To be more specific, we note that in the simplest nonlinear context, a residual network with variable dimensions analog to (2.3) takes the form (see [START_REF] He | Deep residual learning for image recognition[END_REF])

   x k+1 i = Π k x k i + σ w k x k i + b k for k ∈ {0, . . . , N layers -1} x 0 i = x i . (6.1)
Here, contrary to (2.3), we have are given positive integers, called widths of the layers k. One imposes d 0 = d, and Π k ∈ R d k+1 ×d k is a projection/embedding operator which serves to match dimensions. Much like in the fixed width case, we may also write the variable-width ResNet when f is parametrized as in (2.6) or otherwise. The continuous space-time network. It is not immediately obvious how one can see (6.1) as a numerical scheme for some continuous-time dynamical system in the flavor of (2.4). Nevertheless, this can be achieved by viewing the changing dimension over timesteps as an additional (spatial) variable, thus yielding an integro-differential equation in the continuum. To be more precise, for any i ∈ [N ] we consider the scalar integro-differential equation

w k ∈ R d k+1 ×d k and b k ∈ R d k+1 ,
     ∂ t x i (t, x) = σ Ω w(t, x, ξ)x i (t, ξ) dξ + b(t, x) for (t, x) ∈ (0, T ) × Ω x i (0, x) = x in i (x) for x ∈ Ω. (6.2)
Here Ω ⊂ R d Ω is a bounded domain, where d Ω 1. We emphasize that x i (t, x) ∈ R for (t, x) ∈ (0, T ) × Ω, and similarly,

w(t, x, ξ) ∈ R and b(t, x) ∈ R for (x, ξ) ∈ Ω × Ω. The initial datum x in i ∈ C 0 (Ω) is such that there exist {x j } d j=1 ⊂ Ω such that x in i (x j ) = ( x i ) j .
Such a datum can always be found (e.g. by interpolation). The continuum model (6.2) is proposed in [START_REF] Liu | Selection dynamics for deep neural networks[END_REF] where well-posedness is established, and is also suggested in [START_REF] Weinan | A proposal on machine learning via dynamical systems[END_REF] albeit in a slightly different context. We distinguish two typical cases for choosing the shape of Ω as well as d Ω .

• Variable-width ResNets. If in the discretized level, we seek to simply obtain a variable-width residual network such as (6.1) (or even the standard ResNet analog (2.3)), it suffices to consider Ω = (0, 1), thus d Ω = 1. We give more detail on possible possible discretizations in Section 6.1 and Remark 10. • Convolutional Neural Networks. The situation is slightly more delicate in the case of CNNs 9 , which are typically used in tasks arising in computer vision. We provide a proposal covering the continuous-time analog of CNNs with partial generality. Assume that the dataset { x i } i∈[N ] consists of N images: x i ∈ R d 1 ×d 2 ×d ch for any i; here d 1 (resp. d 2 ) denote the number of horizontal (resp. vertical) pixels in the image x i , whereas d ch denotes the number of channels, i.e. the color format (e.g. d ch = 1 for grayscale, d ch = 3 for RGB). In this case, we consider Ω := Ω img × (0, 1), where Ω img ⊂ R 2 is a rectangle. Thus d Ω = 3. Moreover, we assume that the weights w in (6.2) are compactly supported and of a specific convolutional form (as indicated in most works, this is more so a cross-correlation form), namely, for any i, the equation takes the form

∂ t x i (t, x, ζ) = σ 1 0 Ω img w(t, x + ξ, ω, ζ)x i (t, ξ, ω) dξ dω + b(t, x, ζ)
for (t, x, ζ) ∈ (0, T ) × Ω img × (0, 1). We note that the variable x ∈ Ω img denotes a pixel, whereas ζ ∈ (0, 1) is a continuous variable indicating, when discretized,

9
The mathematical theory of structural properties of CNNs in feed-forward form (without skipconnections) is well-established -see for instance [START_REF] Mallat | Understanding deep convolutional networks[END_REF] and the references therein.

the number of extracted features (namely the number of filters). The bias parameter b can be omitted in this case, if desired. One possible way to discretize the above continuous-time model and obtain a CNN-ResNet as in [START_REF] He | Deep residual learning for image recognition[END_REF] is to follow the arguments in Section 6.1, where one would use a time-dependent grid for discretizing with respect to the variable ζ ∈ (0, 1) as well, as the number of filters commonly varies over layers in CNNs. By discretizing Ω img with a "shrinking" or "expanding" timedependent rectangular grid, some effects of padding or pooling (but not maxpooling a priori) may also be considered. However, a full CNN-applicable theory is beyond the scope of this work.

Remark 10. Observe that the continuous space-time model (6.2) (resp. (6.3)) is more general and englobes (2.4) -(2.5) (resp. (2.4) -(2.6)), where only the time variable is considered to be continuous. Indeed, fix d different points {x 1 , . . . , x d } ∈ Ω, and let δ x j denote the Dirac mass centered at x j . For any i ∈ [N ], we consider the initial datum

x in i (x) := d j=1 ( x i ) j δ x j (x) for x ∈ Ω.
We write the weight w as

w(t, x, ζ) := d j=1 d =1 w j, (t)δ x j (x)δ x (ζ) for (t, x, ζ) ∈ (0, T ) × Ω × Ω,
yielding the matrix [w j, (t)] 1 j, d of weights at time t, whereas the bias b(t, x) is written as

b(t, x) := d j=1 b j (t)δ x j (x) for (t, x) ∈ (0, T ) × Ω,
yielding the vector [b j (t)] 1 j d of biases at time t. As x in i , w and b are all linear combinations of Dirac masses, by plugging them in (6.2), we rewrite the integrals as sums, and setting, for any i ∈ [N ],

(x i ) j (t) := Ω x i (t, x) dδ x j (x) for j ∈ [d], we see that (x i ) j solves        ( ẋi ) j (t) = σ d =1 w j, (t) (x i ) (t) + b j (t)
for t ∈ (0, T )

(x i ) j (0) = ( x i ) j .
This is just the j-th equation of the (2.4) -(2.5)

for i ∈ [N ].
Correspondingly for i ∈ [N ] we may consider

     ∂ t x i (t, x) = 1 0 w(t, x, ξ)σ(x i (t, ξ)) dξ + b(t, x) in (0, T ) × Ω x i (0, x) = x in i (x)
in Ω.

(6.3) All of the above discussions also apply for this system.

6.1. From continuous to discrete. The passage from (6.2) to a discrete-time scheme such as (6.1) is not immediately obvious, and to our knowledge has not been presented in the literature. To proceed, it is important to observe the inherent link between the layer k and the width d k in (6.1). This motivates discretizing (6.2) in the spatial variable x ∈ (0, 1) by using a time-dependent grid, which has a different number of nodes d k at each time-step. We give more detail on this in what follows.

Let us demonstrate that (6.2) which reads 10 (we omit the dependence on i for notational simplicity)

     ∂ t x(t, x) = σ 1 0 w(t, x, ξ)x(t, ξ) dξ + b(t, x)
in (0, T ) × (0, 1)

x(0, x) = x in (x) in (0, 1),
where

x in ∈ C 0 ([0, 1]) is such that x in (x j ) = x ,j for some {x j } d j=1 ⊂ [0, 1], can be discretized to read exactly as    x k+1 = Π k x k + σ w k x k + b k for k ∈ {0, . . . , N layers -1} x 0 = x.
(6.4)

Here x k ∈ R d k , w k ∈ R d k+1 ×d k and b k ∈ R d k+1 , with d 0 := d and {d k } N layers k=1
given positive integers, and

Π k ∈ R d k+1 ×d k .
The derivation below is purely for illustrative purposes -an adaptive solver ought to perform better than an adaptation of an Euler scheme as (6.4). Moreover, the subsequent arguments will of course also apply to (6.3). Let t 0 , . . . , t N layers , with t 0 := 0 and t N layers := T, be a given, non-decreasing sequence of time-steps. For simplicity of presentation, let us assume that the time-steps are uniform, namely t k = k∆t with ∆t = T N layers , but more general time-adaptive sequences can be considered. For any k ∈ {0, . . . , N layers }, let us assume that we are given a grid

x j t k d k j=1 ⊂ [0, 1]
which is ordered and uniformly distributed. For simplicity of presentation, in our discussion we will assume that x 1 (t k ) = 0 and x d k (t k ) = 1 for any k. However by means of a time-step-dependent dilation, this restriction may be removed. We note that, not only there might be no overlap of grid nodes over different time-steps, but moreover, the number of grid nodes changes at each time-step k.

We will seek for an appropriate discretization of

∂ t x(t k+1 , x j (t k+1 )) = σ 1 0 w t k+1 , x j (t k+1 ), ξ x(t k , ξ) dξ + b t k+1 , x j (t k+1 )
(6.5) for k ∈ {0, . . . , N layers -1} and j ∈ {1, . . . , d k+1 }. Hence, in view of the preceding discussion, some kind of interpolation may be needed to justify a backward Euler discretization of the time derivative ∂ t x appearing in (6.5) at the grid nodes. 10 The choice of the spatial interval [0, 1] is completely arbitrary -one may of course consider any bounded interval of R.

For any given k ∈ {0, . . . , N layers -1} and j ∈ {1, . . . , d k }, we shall henceforth denote

x k j := x j (t k ),
x k j := x(t k , x k j ). Following through the above discussion, the main issue in writing down a forward difference discretization to ∂ t x(t k+1 , x j (t k+1 )) appears whenever for a given k one has d k = d k+1 , as it is a priori not possible to make sense of the expression x(t k+1 , x j (t k+1 ))-

x(t k , x j (t k )) for j = 1. Indeed, all ι ∈ {2, . . . , d k } are such that x ι (t k ) / ∈ {x j (t k+1 )} d k+1
j=1 , due to the uniformity of the grid.

Let us give an elementary argument for addressing this issue. Given k and given any j ∈ {1, . . . , d k+1 }, there clearly exists ι ∈ {2, . . . ,

d k } such that x k+1 j ∈ [x k ι-1 , x k ι ].
For such indices, we may thus define the linear interpolant

x k j := x k ι + x k ι -x k ι-1 x k ι -x k ι-1 x k+1 j -x k ι . (6.6)
This is nothing but an approximation of the first order Taylor expansion of x(t k+1 , x j (t k+1 )) with respect to the second variable. Using this interpolant, we may consider the simple forward difference

∂ t x(t k+1 , x j (t k+1 )) ≈ x k+1 j -x k j ∆t (6.7)
for any k ∈ {0, . . . , N layers -1} and any j ∈ {1, . . . , d k+1 }. We may now use any Newton-Cotes formula to discretize the integral term in (6.5): for j ∈ {1, . . . , d k+1 }, we write

1 0 w t k+1 , x j (t k+1 ), ξ x(t k , ξ) dξ ≈ d k ι=1 α ι w t k+1 , x j (t k+1 ), x ι (t k ) x(t k , x ι (t k )).
(6.8) Here, α ι > 0 are the corresponding weights of the chosen Newton-Cotes formula.

Let us now define

x k :=    x(t k , x 1 (t k )) . . . x(t k , x d k (t k ))    ∈ R d k , b k :=    b(t k+1 , x 1 (t k+1 )) . . . b(t k+1 , x d k+1 (t k+1 ))    ∈ R d k+1 and w k := α ι w(t k+1 , x j (t k+1 ), x ι (t k )) 1 j d k+1 1 ι d k ∈ R d k+1 ×d k .
The above definitions, as well as (6.7) and (6.8) applied to (6.5), lead us to (6.4), where ∆t has been "omitted" as a factor of the nonlinearity. In view of (6.6), the operator

Π k ∈ R d k+1 × R d k takes the explicit form Π k = d k+1 j=1 1 + x k+1 j -x k ι(j) x k ι(j) -x k ι(j)-1 e j e ι(j) - x k+1 j -x k ι(j) x k ι(j) -x k ι(j)-1 e j e ι(j)-1 , where ι(j) ∈ {2, . . . , d k } is such that x k+1 j ∈ [x k ι(j)-1 , x k ι(j) ], while {e j } d k+1 j=1 and {e j } d k j=1
denote the canonical bases of R d k+1 and R d k respectively. We notice that the matrix Π k only has 2 non-zero elements at every row j ∈ {1, . . . , d k+1 }. This concludes our derivation.

Remark 11 (Generating moving grids). Whilst we have assumed a very simple given time-dependent grid, one may certainly generate more sophisticated moving grids -we refer to [Budd et al., 2009, Section 3] for a comprehensive overview on the existing methods, which have found extensive use in the discretization of partial differential equations manifesting shock waves and/or free boundaries.

6.2. The supervised learning problem. Given a training dataset

{ x i , y i } i∈[N ] with x i ∈ R d for any i ∈ [N ],
just as in the finite dimensional context, we begin by writing the equation satisfied by the stacked vector of states x := [x 1 , . . . , x N ] corresponding to the stacked vector of data

x in := [x in 1 , . . . , x in N ] ,
where each x i is the solution to either (6.2) or (6.3) corresponding to the datum x in i , and control parameters [w, b] which are the same for all i. The stacked continuous space-time neural networks we consider are thus either

     ∂ t x(t, x) = σ Ω w(t, x, ξ)x(t, ξ) dξ + b(t, x) in (0, T ) × Ω x(0, x) = x in (x) in Ω (6.9) or    ∂ t x(t, x) = Ω w(t, x, ξ)σ(x(t, ξ)) dξ + b(t, x) in (0, T ) × Ω x(0, x) = x in (x) in Ω.
(6.10)

Just as in the finite-dimensional case, it is important to note how [w(t, x, ξ), b(t, x)] for (t, x, ξ) ∈ (0, T ) × Ω × Ω enter the systems:

w(t, x, ξ) :=    w(t, x, ξ) . . . w(t, x, ξ)    ∈ R N ×N , b(t, x) :=    b(t, x) . . . b(t, x)    ∈ R N .
6.2.1. Empirical risk minimization. As before, we first consider the regularized empirical risk minimization problem inf [w,b]∈H k (0,T ;U) subject to (6.9) (resp.(6.10))

E(x(T )) + λ [w, b] 2 H k (0,T ;U) , (6.11)
where α > 0 is fixed, k = 0 for (6.10) and k = 1 for (6.9), U := L 2 (Ω × Ω) × L 2 (Ω), and we define the training error as

E(x(T )) := 1 N N i=1 loss(P x i (T ), y i ), (6.12)
where loss ∈ C 0 (R m × Y; R + ) and P : L 2 (Ω) → R m is given. The optimization problem (6.11) admits a solution by the direct method in the calculus of variations.

In view of the rather universal nature of the proof to Theorem 3.1 and Theorem 4.1 in the finite-dimensional case, one may in fact roughly repeat the exact same proofs at most points, replacing throughout the finite dimensional euclidean spaces R dx and R du , by L 2 (Ω) N and U respectively. Whence, we state the infinite-dimensional (partial) analog to Theorem 3.1.

Theorem 6.1. Let λ > 0 be fixed, and let x in ∈ C 0 (Ω) N be such that x in i (x j ) = ( x i ) j . Suppose that P : L 2 (Ω) → R m is any non-zero affine map, and suppose that loss ∈ C 0 (R m × Y; R + ) is such that Assumption 3 is satisfied. Assume that (6.9) (resp. (6.10) with σ positively homogeneous of degree 1) interpolates the dataset {x in i , y i } i∈[N ] in time 1 in the sense of Definition 3.1. For any T 1, let x T ∈ C 0 [0, T ]; L 2 (Ω) N be the unique solution to (6.9) (resp. (6.10)), associated to any global minimizer [w T , b T ] ∈ H k (0, T ; U) of the functional in (6.11), where k = 0 in the case of (6.10) and k = 1 in the case of (6.9). The following properties then hold.

(i) There exists a constant

C = C { x i , y i } i∈[N ] , λ > 0 independent of T such that E(x T (T )) C T .
(ii) There exists a sequence

{T n } ∞ n=1 , with T n > 0 and T n ----→ n-→∞
∞, and some

x • ∈ L 2 (Ω) N with E(x • ) = 0 such that, along a subsequence, E (x Tn (T n )) ----→ n-→∞ 0 and x Tn (T n ) ---- n-→∞ x • weakly in L 2 (Ω) N .
For the sake of completeness, we give a sketch of the proof -by indicating the only changes with respect to that of Theorem 3.1.

Proof of Theorem 6.1. We note that the infinite-dimensional analog of Lemma 7.1 may easily be shown to hold, and one may readily repeat precisely the same arguments as in the proof of Theorem 3.1, replacing R du and R dx by U and L 2 (Ω) N respectively throughout. The only difference occurs in regarding the arguments on strong L 2 -convergence of the sequence of controls in the case k = 1 -in the infinite dimensional case, we may exhibit the Aubin-Lions compactness lemma instead of Rellich-Kondrachov to conclude.

6.2.2. Augmented empirical risk minimization. We similarly consider the augmented supervised learning problem inf [w,b]∈L 2 (0,T ;U) subject to (6.10)

E(x(T )) + 1 N T 0 x(t) -x 2 L 2 (Ω) dt + λ [w, b] 2 L 2 (0,T ;U) , (6.13)
where E is as in (6.12) and loss(•, •) satisfies Assumption 4. We solely concentrate on System (6.10) to avoid a possibly abundance of technical details in the context of BV t L 2 x analysis. Again, P : L 2 (Ω) -→ R m is a surjective map, and

x i ∈ P -1 ({ y i }) ⊂ L 2 (Ω) for i ∈ [N ]
are arbitrary, but fixed. As expected, the analog exponential decay result holds for (6.13).

Theorem 6.2 (Exponential stability). Fix λ > 0, let P ∈ Lip(L 2 (Ω); R m ) be any given surjective map, and let x ∈ L 2 (Ω) N with x i ∈ P -1 ({ y i }) for i ∈ [N ] be arbitrary but fixed. Suppose that (6.10) is controllable with linear cost in some time T 0 > 0 in the sense of Definition 4.1. Then, there exists T * > 0 and constants

C = C { x i , y i } i∈[N ] , λ, N, α, P > 0 and µ = µ { x i , y i } i∈[N ]
, λ, N, α > 0 such that for any T T * , any pair of parameters [w T , b T ] ∈ L 2 (0, T ; U) solving the minimization problem (6.13), and the corresponding unique solution x T (•) to (6.10) satisfy

w T (t) + b T (t) C for a.e. t ∈ [0, T ] and E(x T (t)) + x T (t) -x L 2 (Ω) C e -µt for all t ∈ [0, T ].
The proof is omitted and left to the reader, as it follows precisely the same arguments as that of Theorem 4.1.

7. Proofs 7.1. Proof of Theorem 3.1. We note that both (3.2) and (3.3) can be written in the compact form

ẋ(t) = f([w(t), b(t)], x(t)) in (0, T ) x(0) = x 0 ∈ R dx , (7.1) with f([0, 0], x) = 0, f([αw, αb], x) = αf([w, b], x) for α > 0. (7.
2) We will refer to u := [w, b] as the control of the ODE system, in accordance with control theory vocabulary. We begin with following short but key lemma.

Lemma 7.1. Let T 0 > 0 and [w T 0 , b T 0 ] ∈ L 1 (0, T 0 ; R du ) be given, and let x T 0 (•) be the unique solution to

ẋT 0 (t) = f ([w T 0 (t), b T 0 (t)] , x T 0 (t)) in (0, T 0 ) x T 0 (0) = x 0 ∈ R dx , (7.3)
(i.e. (7.1) on (0, T 0 )) with f as in either (3.3) or (3.2), thus satisfying (7.2). Let T > 0, and define

w T (t) := T 0 T w T 0 t T 0 T , b T (t) := T 0 T b T 0 t T 0 T for t ∈ [0, T ],
and

x T (t) := x T 0 t T 0 T for t ∈ [0, T ].
Then x T (•) is the unique solution to (7.1) (with the same f as in

(7.3)) associated to [w T , b T ].
This sort of time-scaling in the context of driftless control affine systems is commonly used in control theoretical contexts -a canonical example is the proof of the Chow-Rashevskii controllability theorem, see [Coron, 2007, Chapter 3, Section 3.3]. We sketch the short proof for completeness.

Proof of Lemma 7.1. Since x T 0 is the solution to (7.3), the change of variable τ = s T T 0 as well as (7.2), we have

x T (t) := x T 0 t T 0 T = x 0 + t T 0 T 0 f ([w T 0 (s), b T 0 (s)] , x T 0 (s)) ds = x 0 + t 0 T 0 T f w T 0 τ T 0 T , b T 0 τ T 0 T , x T 0 τ T 0 T dτ = x 0 + t 0 f ([w T (τ ), b T (τ )] , x T (τ )) dτ.
It follows that x T solves (7.1), and we conclude by uniqueness.

We will also need the following lemma.

Lemma 7.2 (Compactness of the flow). Let T > 0 be fixed. The maps (i)

Φ T : [w, b] -→ x(•) mapping L 2 (0, T ; R du ) to C 0 ([0, T ]; R dx ) where x(•) solves (3.3), (ii) Φ T : [w, b] -→ x(•) mapping L 2 (0, T ; R du ) ∩ BV(0, T ; R du ) to C 0 ([0, T ]; R dx ) where x(•) solves (3.2), (iii) Φ T : [w, b] -→ x(•) mapping H 1 (0, T ; R du ) to C 0 ([0, T ]; R dx ) where x(•) solves (3.2), are all compact.
We postpone the proof to the appendix. We are now in a position to prove the main result.

Proof of Theorem 3.1. We will henceforth, for notational convenience, extensively make use of the notation u := [w, b]. We will focus on the neural ODE (3.3) and hence k = 0. The case (3.2) and k = 1 follows exactly the same arguments, and we will comment on the key differences at the end of the proof. Part 1. We begin by showing E (x T (T )) T -1 (7.4) uniformly in T . By the interpolation assumption, there exists some u 1 ∈ L 2 (0, 1; R du ) such that the associated solution x 1 to (3.3) on [0, 1] satisfies E(x 1 (1)) = 0. Using the optimality of u T and the scaling relations from Lemma 7.1, we obtain

J λ,T (u T ) = E (x T (T )) + λ u T 2 L 2 (0,T ;R du ) E(x 1 (1)) + λ T u 1 2 L 2 (0,1;R du )
for all T > 0. Since E(x 1 (1)) = 0 by the interpolation assumption, the above inequality implies

0 E (x T (T )) λ T u 1 2 L 2 (0,1;R du ) (7.5)
for all T > 0. Estimate (7.5) clearly implies (7.4). Part 2. We now look to prove (3.6). To this end, we will look to show that {x T (T )} T >0 is a bounded subset of R dx . This will allow us to extract a converging sequence, whose limit will be shown to lie in {E = 0}.

For any T > 0, set

u aux (t) := 1 T u 1 t T for t ∈ [0, T ].
We argue similarly as in Part 1. Making use of Lemma 7.1 once again, and since E(x 1 (1)) = 0, we see that

J λ,T (u aux ) = E x 1 (1) + λ T u 1 2 L 2 (0,1;R du ) = λ T u 1 2 L 2 (0,1;R du ) . (7.6)
Using the optimality of u T , one sees that

J λ,T (u aux ) J λ,T (u T ) λ u T 2 L 2 (0,T ;R du ) . (7.7) that J λ,Tn u † n = λ T n u † 2 L 2 (0,1;R du ) < E (x Tn (T n )) + λ T n u n 2 L 2 (0,1;R du ) = J Tn (u Tn ) ,
which contradicts the fact that u Tn minimizes J Tn . Hence, (7.11) holds, and {u n } ∞ n=1 is bounded in L 2 (0, 1; R du ). Consequently, by the Banach-Alaoglu theorem, there exists

u * = [w * , b * ] ∈ L 2 (0, 1; R du ) such that u n u * weakly in L 2 (0, 1; R du ),
along some subsequence as n -→ ∞. Moreover, using the properties of equation (3.3) (Lemma 7.2), we deduce that the trajectory x n associated to u n satisfies

x n -→ x * strongly in C 0 ([0, 1]; R dx ) (7.12)

as n -→ ∞, where x * is the solution to (3.3) on [0, 1], associated to u * . On another hand, note that by Lemma 7.1, x Tn (t) = x n ( t Tn ) for t ∈ [0, T n ], whence x Tn (T n ) = x n (1) and thus, combining (7.12) and (7.9), we see that x * (1) = x • . Consequently, u * is a control such that E(x * (1)) = E(x • ) = 0, thus satisfying the constraint in (7.10). In view of this, we may also use (7.11) and the weak lower semicontinuity of the L 2 -norm to write

u † L 2 (0,1;R du ) u * L 2 (0,1;R du ) lim inf n-→∞ u n L 2 (0,1;R du ) lim n-→∞ u n L 2 (0,1;R du ) lim sup n-→∞ u n L 2 (0,1;R du ) u † L 2 (0,1;R du ) , (7.13) clearly implying that lim n-→∞ u n L 2 (0,1;R du ) = u * L 2 (0,1;R du ) .
Hence, as weak convergence and convergence of the norms in L 2 implies strong convergence in L 2 , we deduce that

u n -→ u * strongly in L 2 (0, 1; R du )
along some subsequence as n -→ ∞. Moreover, from (7.13) we deduce that, since u † is a solution to (7.10) and since u * satisfies the constraints therein, u * is a solution to (7.10) as well, which concludes the proof for (3.3) and k = 0.

In the case (3.2) and k = 1, one may clearly repeat the above reasoning, replacing L 2 (0, T ; R du ) by H 1 (0, T ; R du ) throughout, with some key additions.

In Part 1, we first note that instead of (7.6), one has

J λ,T (u aux ) = E(x 1 (1)) + λ T u 1 2 L 2 (0,1;R du ) + λ T 3 u1 2 L 2 (0,1;R du ) = λ T u T 0 2 L 2 (0,1;R du ) + λ T 3 u1 2 L 2 (0,1;R du ) .
This is not an impediment to (7.7), which remains true, and one can clearly deduce that {x T (T )} T >0 is bounded as well. Similarly, (7.5) holds with a bound of the form

0 E(x T (T )) λ T u 1 2 L 2 (0,1;R du ) + λ T 3 u1 2 L 2 (0,1;R du ) .
Whence the remainder of parts 1 and 2 hold in this context as well.

In Part 3, we emphasize the sole key difference between (3.3) and (3.2) -the weak L 2convergence of {u n } ∞ n=1 is a priori not sufficient to entail the strong convergence in (7.12) in the case of (3.2). However, by the Rellich-Kondrachov compactness theorem, the weak H 1 -convergence of {u n } ∞ n=1 implies a strong L 2 -convergence along a subsequence, which would yield (7.12) by arguing just as in the proof of Lemma 7.2.

This concludes the proof.

7.2. Proof of Theorem 3.2. The proof closely follows the lines of that just above. Let us consider k = 1, since the case k = 0 is equivalent to Theorem 3.2. We present minimal details for completeness.

Proof of Theorem 3.2. We again make use of the notation u := [w, b]. We first show

E(x λ (T )) λ (7.14)
uniformly in λ > 0 -we argue as in the proof of Theorem 3.1 just above, exhibiting, by the interpolation assumption, parameters u 1 ∈ L 2 (0, 1; R du ) such that E(x 1 (1)) = 0. We may obtain an estimate like (7.5) and conclude. Now, the same arguments as in Part 2 of the proof of Theorem 3.1 may be used to deduce that {x λ (T )} λ>0 is a bounded subset of R d , and hence there exists a sequence {λ n } ∞ n=1 of positive numbers with λ n 0 as n -→ ∞ and some x • ∈ R dx such that

x λn (T ) ----→ n-→∞ x • .
Using (7.14) we deduce that E(x • ) = 0. Finally, the proof of the last fact is identical to that done for Theorem 3.1, so we omit it. 7.3. Proof of Theorem 3.3. We now provide a proof of our main result in the context of classification tasks.

Proof of Theorem 3.3. We recall that, by assumption,

x 0 i := Q x i 0 for i ∈ [N ].
Now let w, b ∈ H k (0, T 0 ; R du ) be a pair of parameters which separates the training dataset x 0 i , y i i∈[N ] with respect to P in time T 0 > 0, i.e., such that the solution x = [ x 1 , . . . , x N ] to (3.2), with initial condition x 0 = [x 0 1 , . . . , x 0 N ], corresponding to w, b , satisfies

min i∈[N ]      P x i (T 0 ) y i -max j∈[N ] j = y i P x i (T 0 ) j      =: γ > 0. (7.15)
Now fix an arbitrary α ∈ 0, 1 2 , and, for any T > 0, define

w † (t), b † (t) :=        2T 0 T w t 2T 0 T , b t 2T 0 T for t ∈ 0, T 2 T α-1 [Id d , 0 d ] for t ∈ T 2 , T ,
where Id d is the identity matrix in R d×d and 0 d is the zero vector in R d . By virtue of the scaling in Lemma 7.1, for t ∈ T 2 , T , the trajectories x † = x † 1 , . . . , x † N associated to w † , b † are given by the solution to

       ẋ † i (t) = σ T α-1 x † i (t) for t ∈ T 2 , T x † i T 2 = x i (T 0 ) . (7.16)
Moreover, since σ(x) = max{x, 0} 0, the right hand side in (3.2) is nonnegative. Using the assumption that the initial conditions are of the form x 0 i = Q x i 0, it follows that x i (T 0 ) 0 for all i ∈ [N ]. We can therefore drop σ in (7.17) and deduce that P

x † i (t) solves        d dt P x † i (t) = T α-1 P x † i (t) for t ∈ T 2 , T P x † i T 2 = P x i (T 0 ) .
(7.17)

Hence, we have

P x † i (t) = P x i (T 0 )e T α-1 (t-T /2) , for all t ∈ T 2 , T .
Now, using the definition of the cross-entropy loss and the margin γ in (7.15), we compute, for any i ∈ Then, we can estimate

E x † (T ) log 1 + (m -1) exp -γ exp T α 2 . (7.18)
On the other hand, using the definition of w † , b † , we deduce

w † , b † 2 H 1 (0,T ;R du ) = w † , b † 2 H 1 (0, T 2 ;R du ) + w † , b † 2 H 1 ( T 2 ,T ) C 1 T + C 2 T 2(α-1) T,
for some constants C 1 , C 2 > 0 depending only on λ, T 0 and w, b . From this estimate, together with (7.18), we obtain, for T > T 0 ,

J λ,T w † , b † log 1 + (m -1) exp -γ exp T α 2 + CT 2α-1 ,
for some constant C > 0 depending on λ, T 0 , w, b , but independent of T . Using the above estimate, we may conclude from the optimality of [w T , b T ], as

E (x T (T )) J λ,T (w T , b T ) J λ,T w † , b † log 1 + (m -1) exp -γ exp T α 2 + CT 2α-1 .
7.4. Proof of Theorem 4.1. The proof of Theorem 4.1 is rather involved -we will first require the following lemmas. We will focus on System (3.2), but with small adaptations the proofs (except the exponential stability estimate for the parameters) can be shown for systems of the form (3.9). System (3.3) is addressed in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF]].

Lemma 7.3 (Grönwall-like estimate). Let T > 0 be given, and let x ∈ R dx . For any [w, b] ∈ L 1 (0, T ; R du ) and x 0 ∈ R dx , let x ∈ C 0 ([0, T ]; R dx ) denote the unique solution to (3.2), noting (3.1). The following facts then hold.

(i) There exist a couple of constants

C 1 = C 1 σ, x , max 1, x 0 -x , N > 0 and C 2 = C 2 (σ, x , N ) > 0 independent of T such that x(t) -x C x 0 -x + [w, b] L 1 (0,T ;R du ) + x -x L 2 (0,T ;R dx )
holds for all t ∈ [0, T ], where

C := C 1 max 1, [w, b] L 1 (0,T ;R du ) exp C 2 w L 1 (0,T ;R d×d ) . (ii) If moreover [w, b] ∈ L 2 (0, T ; R du ), then x(t) -x C x 0 -x + [w, b] L 2 (0,T ;R du ) + x -x L 2 (0,T ;R dx )
also holds for all t ∈ [0, T ], where

C := C 1 max 1, [w, b] L 2 (0,T ;R du ) exp C 2 w L 2 (0,T ;R d×d ) .
The proof of Lemma 7.3 may be found in the Appendix. Now suppose x τ• ∈ R dx is given. Let T > 0 and τ • ∈ [0, T ) be fixed, and consider the cost functional

J τ•,T (w, b) := 1 N T τ• x(t) -x 2 dt + λ [w, b] 2 L 2 (τ•,T ;R du ) , (7.19)
with x(•) being the solution to (7.20) with f as the dynamics in (3.3) or (3.2). We then have the following result.

ẋ(t) = f ([w(t), b(t)], x(t)) in (τ • , T ) x(τ • ) = x τ• ,
Lemma 7.4 (Uniform estimate of optimal trajectories). Fix λ > 0, let P ∈ Lip(R d ; R m ) be any given surjective map, and let x ∈ R dx with x i ∈ P -1 ({ y i }) for i ∈ [N ] be arbitrary but fixed. Let r > 0 be fixed. Suppose that (3.3) (resp. (3.2) with σ 1-homogeneous) is controllable with linear cost in the sense of Definition 4.1. Then, there exists a constant C = C (σ, x , λ, N, d, r) > 0 such that for all T > 0, τ • ∈ [0, T ) and x τ• ∈ R dx such that x τ•x r, any pair of parameters [w T , b T ] ∈ L 2 (τ • , T ; R du ) minimizing J τ•,T and corresponding solution x T (•) to (7.20), with f as in (3.3) (resp (3.2) with σ 1-homogeneous), are such that

[w T , b T ] L 2 (τ•,T ;R du ) + x T -x L 2 (τ•,T ;R dx ) + x T (t) -x C x τ• -x holds for all t ∈ [τ • , T ].
The proof of Lemma 7.4 follows precisely the arguments given in the first step of the proof of Theorem 4.2, in which one solely needs to change the initial time (0 by τ • ) and datum (x 0 by x τ• ), the BV-norms by L 2 , and note the precise form of the constant C in (7.30). We give a brief sketch in the appendix for the sake of clarity. We will also make use of the following short observation.

Lemma 7.5. Let T > 0 and let [w T , b T ] ∈ L 2 (0, T ; R du ) be a pair of minimizers to J T defined in (4.2), and denote by x T the corresponding solution to (7.20) on [0, T ] with

x T (0) = x 0 . Let τ • ∈ [0, T ) be given. Then [w T , b T ]| [τ•,T ] minimize J τ•,T defined in (7.19) for System (7.20) with initial data x τ• = x T (τ • ).
We refer to the appendix for a proof.

Remark 12. We note that our proof is specific to L 2 -regularization and does not transfer to BV-regularization a priori, due to the fact that the BV norm may see the singularities in discontinuous parameters. Therein lies the main impediment to ensuring the validity of Theorem 4.1 for BV-regularized problems by means of our strategy.

Before concluding this section with a proof of Theorem 4.1, we also state the following key lemma.

Lemma 7.6. Let X be a Banach space, T > a 0 and f ∈ C 0 ([a, T ]; X). For any τ Ta, there exists

t 1 ∈ [a, a + τ ) such that f (t 1 ) X f L 2 (a,T ;X) √ τ .
The proof may be found in the appendix. We may now provide the proof of Theorem 4.1.

Proof of Theorem 4.1. We shall concentrate on the neural ODE (3.2). The proof of the full result for (3.3) is identical to that presented in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF]. We split the proof in two parts.

Part 1: Stability estimates for E(x T (t)) + x T (t)x . For r := x 0x , denote by C 1 = C 1 (σ, x , λ, N, d, r) > 0 the universal constant given by Lemma 7.4. Let τ > max C 4 1 , C 2 1 be fixed and let T τ + 1. First, note that by Lemma 7.4 (with τ • = 0, r := x 0x and x τ• = x 0 ), we have

x T -x L 2 (0,T ;R dx ) + x T (t) -x C 1 x 0 -x for t ∈ [0, T ]. (7.21)
Now note that for t ∈ [0, τ + 1], the desired exponential stability estimate can easily be obtained since the length of the time interval is independent of T . Indeed, from (7.21) we get

x T (t)x C 1 x 0x exp(t) exp(-t)

C 1 x 0x exp(τ + 1) exp(-t),

and, since x i ∈ P -1 ({ y i }) for i ∈ [N ], using this estimate we also find

E(x T (t)) c N N i=1 P x T,i (t) -y i α c N P α N i=1 x T,i (t) -x i α C α 1 x 0 -x α exp(α(τ + 1)) exp(-αt).
Thus, only the case t ∈ [τ + 1, T ] remains. To do so, we proceed in three steps.

Step 1). Preparation. Since τ T , using Lemma 7.6 and then Lemma 7.4 (with τ • = 0, r := x 0x and x τ• = x 0 ) we see that there exists τ Step 2). Bootstrap. We iterate (7.23) and show that for any n ∈ N satisfying n T /τ, the following estimate holds:

• ∈ [0, τ ) such that x T (τ • ) -x x T -x L 2 (0,T ;R dx ) √ τ C 1 √ τ x 0 -x . ( 7
x T (t) -x C 1 x T (τ • ) -x C 2 1 √ τ x 0 -x , ( 7
x T (t) -x C 2 1 √ τ n x 0 -x for t ∈ [nτ, T ]. (7.24)
We proceed by induction -the case n = 1 holds by (7.23). Thus suppose that (7.24) holds at some stage n 2 and suppose that n + 1 T /τ. Now the parameters [w T , b T ]| [nτ,T ] minimize J nτ,T by Lemma 7.5, and so we can apply Lemma 7.4 (with τ • = nτ , x τ• = x T (nτ ), and with r = x 0x as C 2 1/ √ τ < 1 and (7.24) is assumed to hold) in combination with Lemma 7.6 (as n + 1 T /τ clearly implies that τ Tnτ ) to deduce that there exists some time t 1 ∈ [nτ, (n + 1)τ ) such that

x T (t 1 ) -x x T -x L 2 (nτ,T ;R dx ) √ τ C 1 √ τ x T (nτ ) -x .
We once again apply (7.24) in the inequality just above to deduce

x T (t 1 ) -x C 1 √ τ C 2 1 √ τ n x 0 -x . (7.25)
We then apply Lemma 7.4 for a second time (with τ • = t 1 , x τ• = x T (t 1 ) and with r = x 0x , as we have (7.25) with max{C1, C 2 1 } / √ τ < 1) and use (7.25) to find

x T (t) -x C 1 x T (t 1 ) -x C 2 1 √ τ n+1 x 0 -x ,
which holds for all t ∈ [t 1 , T ]. Clearly, as t 1 < (n + 1)τ , the above estimate also holds for all t ∈ [(n + 1)τ, T ]. This completes the proof of (7.24).

Step 3). Conclusion. Suppose t ∈ [τ + 1, T ] is arbitrary and fixed. Set n(t) := t τ +1 ; note that n(t) 1, t n(t)τ and n(t) T /τ. We may then apply (7.24) to find that

x T (t) -x C 2 1 √ τ n(t)
x 0x .

Since τ > C 4 1 and n(t) t τ +1 -1, from the above inequality we infer

x T (t) -x exp -n(t) log √ τ C 2 1 x 0 -x √ τ C 2 1 exp   - log √ τ C 2 1 τ + 1 t   x 0 -x .
The desired exponential stability estimate for x T (t)x thus also holds for

t ∈ [τ +1, T ], with µ := log √ τ C 2 1 τ +1 > 0 and C := √ τ C 2 1
x 0x . Note that, arguing as in the previous case,

E(x T (t)) 1 N N i=1 P x T,i (t) -y i α C α exp(-αµt).
This concludes the proof of the first part.

Part 2: Stability estimate for the parameters. The stability estimate for the parameters in the setting of (3.2) closely follows the proof presented in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF] -we give a sketch of the main ideas. We henceforth interchange between the notation u := [w, b] and [w, b] to ease the reading. Fix an arbitrary t ∈ [0, T ) and 0 < h 1, so that t + 2h 2 + 2h ∈ [0, T ], and set

u aux (s) :=                              u T (s) for s ∈ [0, t] 1 2 u T t + s -t 2 for s ∈ (t, t + 2h 2 ] h + 2 2 u T h + 2 2 s - h + 2 2 (t + 2h 2 ) + t + h 2 for s ∈ (t + 2h 2 , t + 2h 2 + 2h] u T (s) for s ∈ (t + 2h 2 + 2h, T ].
From Lemma 7.1 that the solution x aux to (3.2) associated to the above pair is precisely

x aux (s) :=                          x T (s) for s ∈ [0, t] x T t + s -t 2 for s ∈ (t, t + 2h 2 ] x T h + 2 2 s - h + 2 2 (t + 2h 2 ) + t + h 2 for s ∈ (t + 2h 2 , t + 2h 2 + 2h] x T (s) for s ∈ (t + 2h 2 + 2h, T ].
This specific construction is in particular done as to ensure that x aux (T ) = x T (T ) and so E(x aux (T )) = E(x T (T )). Now by several straightforward computations (which may be found in [START_REF] Esteve | Turnpike in Lipschitznonlinear optimal control[END_REF]) we deduce that

J T (u aux ) = E(x T (T )) + 1 N T 0 x T (s) -x 2 ds + 1 N t+h 2 t x T (s) -x 2 ds + 1 N 2 h + 2 -1 t+2h 2 +2h t+h 2 x T (s) -x 2 ds + λ T 0 u T (s) 2 ds - λ 2 t+h t u T (s) 2 ds + λ h 2 t+h 2 +2h t+h 2 u T (s) 2 ds E(x T (T )) + 1 N T 0 x T (s) -x 2 ds + 1 N t+h t x T (s) -x 2 ds + λ T 0 u T (s) 2 ds - λ 2 t+h t u T (s) 2 ds + λ h 2 t+h 2 +2h t+h 2 u T (s) 2 ds.
The above identity combined with the optimality inequality J T (u T ) J T (u aux ) leads us to

λ 2 t+h t u T (s) 2 ds 1 N t+h t x T (s) -x 2 ds + λ h 2 t+h 2 +2h t+h 2 u T (s) 2 ds.
Using the exponential stability estimate for x T (•)x , we may find

λ 2 1 h t+h t u T (s) 2 ds 1 N 1 h t+h t x T (s) -x 2 ds + λ 2 t+h 2 +2h t+h 2 u T (s) 2 ds C N 1 h t+h t e -2µs ds + λ 2 t+h 2 +2h t+h 2 u T (s) 2 ds C N e -2µt + λ 2 t+h 2 +2h t+h 2 u T (s) 2 ds.
By using the Lebesgue dominated convergence theorem (applied to the integrable function s → u T (s) 2 1 (t+h 2 ,t+h 2 +2h) (s)) the second integral in the estimate just above goes to 0 when h 0. Hence, by applying the Lebesgue differentiation theorem in the estimate just above, we deduce that

[w T (t), b T (t)] 2 = u T (t) 2 = lim h 0 1 h t+h t u T (s) 2 ds C 2N λ e -2µt
for a.e. t ∈ [0, T ], as desired.

7.5. Proof of Theorem 4.2.

Proof of Theorem 4.2. We will interchange the notations [w, b] and u := [w, b] for simplicity, and we split the proof in two parts. We henceforth set r := x 0x .

Part 1: Uniform estimates. We shall first establish uniform-in-T estimates -we find some C > 0 independent of T such that

x T (t) -x + x T -x L 2 (0,T ;R dx ) + [w T , b T ] BV(0,T ;R du ) C (7.26)
whenever T 1 and for t ∈ [0, T ]. By virtue of the controllability assumption, there exist parameters u † ∈ C 0 ([0, 1]; R du ) ∩ BV([0, 1]; R du ) such that the corresponding solution x † (•) to (3.2) on (0, 1) satisfies x † (1) = x. By integrating the equation satisfied by x † (•), we see that for t ∈ [0, 1],

x † (t) -x x 0 -x + c(σ) t 0 w † (s) x † (s) -x ds + x t 0 w † (s) ds + t 0 b † (s) ds
where c(σ) > 0 is the Lipschitz constant of σ. By virtue of Grönwall's inequality and (4.5), we deduce that

N -1 /2 x † (t) -x C 1 x 0 -x exp C 1 x 0 -x (7.27)
for t ∈ [0, 1] and for some constant C 1 = C 1 (σ, x , r) > 0 independent of T (as well as and N , and only depends on x 0 via r). We now set (recall that T 1)

u aux (t) := u † (t) in (0, 1) 0 in (1, T ),
and we denote by x aux (•) the corresponding solution to (3.2). One notes that u aux ∈ L 2 (0, T ; R du ) ∩ BV(0, T ; R du ), as the jump at t = 1 only accounts to a Dirac mass. Moreover, x aux (t) = x for t ∈ [1, T ] and thus also P x aux i (T ) = y i for i ∈ [N ], which implies E(x aux (T )) = 0. By using the optimality inequality J T (u T ) J T (u aux ) and the fact that E(x aux (T )) = 0 and E 0, we thus deduce that

N -1 x T -x 2 L 2 (0,T ;R dx ) + λ u T 2 BV(0,T ;R du ) N -1 x aux -x 2 L 2 (0,T ;R dx ) + λ u aux 2 BV(0,T ;R du ) N -1 x † -x 2 L 2 (0,1;R dx ) + 2λ u † 2 BV(0,1;R du ) + 2λ(d 2 + d) d 2 +d j=1 u † j (1) 2 , (7.28)
where we used the vectorized form of u † in the second component of the BV-norm whilst maintaining the same notation (recall that d u := d × (d + 1)). Using (7.27) and (4.5), we find

N -1 x † -x 2 L 2 (0,1;R dx ) +2λ u † 2 BV(0,1;R du ) + 2λ(d 2 + d) d 2 +d j=1 u † j (1) 2 C 2 x 0 -x 2 exp 2C 1 x 0 -x :=C 2 0 (7.29)
for some constant C 2 = C 2 (σ, x , r, d, max{1, λ}) > 0. Combining (7.28) and (7.29), and recalling Lemma 7.3, we may conclude that (7.26) holds with

C := C 3 max 1, C 0 λ exp C 4 C 0 λ x 0 -x (7.30) for some constants C 3 = C 3 (σ, x , r, N, λ) > 0 and C 4 = C 4 (σ, x , N ) > 0.
Part 2: Conclusion. We note that the desired stability estimates thus follow from (7.26), and, since x i ∈ P -1 ({ y i }) for i ∈ [N ], we also have

E(x T (t)) N -1 N i=1 P x T,i (t) -y i α N -1 x T (t) -x α .
We conclude the proof by noting that the convergence of averages follows by dividing both sides in (7.26) by T , using the estimate just above for the training error term, and letting T -→ ∞.

7.6. Proof of Proposition 5.1. The proof of Proposition 5.1 is a straightforward Grönwall argument. We sketch it for completeness.

Proof of Proposition 5.1. For simplicity of presentation but without any loss of generality, we will henceforth concentrate on system (3.3). For any t ∈ [0, T ], i ∈ [N ] and j ∈ [N ], we have

x i (t) -x j (t) = x 0 i -x 0 j + t 0 w(τ ) σ(x i (τ )) -σ(x j (τ )) dτ.
Using the Lipschitz character of σ, we get

x i (t) -x j (t) x 0 i -x 0 j + t 0 w(τ ) σ (x i (τ )) -σ (x j (τ )) dτ x 0 i -x 0 j + C(σ) t 0 w(τ ) x i (τ ) -x j (τ ) dτ.
We apply the Grönwall inequality with the effect of

x i (t) -x j (t) exp C(σ) t 0 w(τ ) dτ x 0 i -x 0 j .
We evaluate the above expression at final time t = T to obtain

x 1 i -x 1 j exp C(σ) T 0 w(τ ) dτ x 0 i -x 0 j , for some x 1 i ∈ P -1 ({ y i }) and x 1 j ∈ P -1 ({ y j }), whence exp C(σ) T 0 w(τ ) dτ x 1 i -x 1 j x 0 i -x 0 j .
By taking the log on both sides we obtain (5.2). 7.7. Proof of Theorem 5.1. The following short functional analysis lemma will be of use in the proof of Theorem 5.1. We omit the proof, which follows by using the open mapping theorem (see e.g. [Brezis, 2010, Theorem 2.6, pp. 35]).

Lemma 7.7. Let H 1 and H 2 be two real Hilbert spaces. Let Λ : H 1 -→ H 2 be a linear, bounded and surjective operator. Then

Γ : H 2 -→ H 1 y -→ arg min x∈Λ -1 ({y})
x 2

H 1 is linear and bounded.

Proof of Theorem 5.1. Inspired by the techniques in [START_REF] Coron | Global steady-state controllability of onedimensional semilinear heat equations[END_REF][START_REF] Pighin | Controllability under positivity constraints of semilinear heat equations[END_REF], we define the continuous arc

γ : [0, 1] -→ R dx s -→ (1 -s)x 0 + sx 1 .
By assumption, σ x 1 1 , . . . , σ x 1 i , . . . , σ x 1 N is a linearly independent system of vectors in R d for any s ∈ [0, 1]. Thus, by using the continuity of γ, there exists an η > 0, such that whenever x 1x 0 η, σ (γ 1 (s)) , . . . , σ (γ i (s)) , . . . , σ (γ N (s)) (7.31) is also a system of linearly independent vectors in R d for any s ∈ [0, 1]. Following the framework of Lemma 7.7, for any s ∈ [0, 1], define

Λ s : R d×d -→ R dx w -→ wσ (γ (s)) .
By the linear independence of the system of vectors (7.31), Λ s is surjective for any s ∈ [0, 1]. Hence, using Lemma 7.7, we see that

Γ s : R dx -→ R d×d y -→ arg min w∈Λ -1 s ({y})
w , is a linear and bounded operator for any s ∈ [0, 1], and, since (7.31) is independent and the arc γ is continuous, {Γ s } s∈[0,1] is uniformly bounded in operator norm:

Γ s L(R dx ;R d×d ) C (7.32) for some C > 0 independent of T > 0. Now, for t ∈ [0, T ], set w(t) := Γ st x 1 -x 0 T , (7.33)
with s t := t T . Note that for any t ∈ [0, T ], the vector w(t) ∈ R d×d solves the linear system of equations w(t)σ (x i (t)) = ẋi (t) for i ∈ [N ], where

x(t) := γ t T = 1 - t T x 0 + t T x 1 . Hence, x(t) solves      ẋi (t) = w(t)σ(x i (t)) for t ∈ (0, T ) x i (0) = γ(0) = x 0 i x i (T ) = γ(1) = x 1 i , for any i ∈ [N ]
. This thus demonstrates the existence of a control w steering the stacked dynamics from x 0 to x 1 in time T .

Let us conclude by showing that w satisfies the stated estimate. By the definition of w in (7.33) as well as (7.32), for any t ∈ [0, T ] we have

w(t) = Γ t x 1 -x 0 T C T x 1 -x 0 ,
as desired.

Remark 13 (Removing the smallness assumption). One could perhaps adapt the argument in the proof of Theorem 5.1 (given just below) to obtain a global result, assuming the existence of a continuous arc γ linking x 0 and x 1 such that σ (γ 1 (s)) , . . . , σ (γ i (s)) , . . . , σ (γ N (s))

is a system of linearly independent vectors in R d for any s ∈ [0, 1]. Problems arise however whenever this condition is not satisfied. In any case, in view of the uniqueness results for ODEs and Proposition 5.1, we have to assume that x 0 i = x 0 j and x 1 i = x 1 j , for i = j.

Concluding remarks

In this work, we have addressed the impact of the final time horizon T in general learning problems for neural ODEs.

• In the empirical risk minimization problem with a Tikhonov (L 2 or H 1 ) parameter regularization, we concluded via Theorem 3.1 -Theorem 3.3 that when T is large enough, the obtained optimal/trained parameters for neural ODEs are such that the corresponding trajectories reach zero training error with a quantitative rate (thus, stipulate an approximation property of the trained model with respect to T ), whilst doing so with the least oscillations possible. In the associated discrete-time, residual neural network setting, this result indicates that adding more layers before training would guarantee the optimal trajectories approach the zero training error regime, but do so without overfitting. In more practical terms, to ensure that the global minimizer is near zero training error, while training, one could systematically increase the time horizon T whilst keeping the regularization parameter λ > 0 fixed. Moreover, roughly the same conclusions hold when T is fixed and λ is rendered small, thereby linking our insights with the literature on the regularization path limit for various machine learning models.

• To obtain better quantitative estimates on the time horizon (and thus, number of layers when the time-step is fixed, e.g. in ResNets) required to be ε-close to the zero training error regime, for a given tolerance ε > 0, we introduced a minimization problem wherein we added a tracking term which regularizes the state trajectories over the entire time horizon. In Theorem 4.1, we show that the training error and the optimal parameters are at most of the order O e -µt for all t ∈ [0, T ]. This result, along with numerical experiments, demonstrates a strong approximation rate of the trained neural ODE flow (which ought to be compared with universal approximation results, in which, a key caveat is that there is no scalable method to compute the theoretically guaranteed parameters), with parameters which are exponentially small, and could thus stipulate that the flow would tend to oscillate little. Moreover, the exponential decay estimate also ensures that T need not be chosen too large to render the training error small. 8.1. Outlook. We present a list of questions and topics which would be complementary to our work.

Generalization bounds. To complement our analytical study on the long time horizon/large layer regime, it would be of interest to provide generalization error bounds for the limiting, least L 2 -norm parameters in the interpolation regime obtained in Theorem 3.1, via, for instance, commonly used metrics such as the VC dimension [START_REF] Vapnik | The nature of statistical learning theory[END_REF] or Rademacher complexity [START_REF] Bartlett | Rademacher and Gaussian complexities: Risk bounds and structural results[END_REF]. Such studies are, to the best of our knowledge, not done in the context of models such as neural ODEs.

Stability estimate for (4.7). We provided a proof of the exponential stability estimate of the training error and optimal parameters in the context of 2 -like losses, and without regularizing the output P x i (t) but rather the features x i (t) over all time/layer t ∈ [0, T ].

We could stipulate that, whenever P is Lipschitz (and possibly real analytic) and such that the training error attains its minimum (e.g. when P is a matrix, or a matrix composed with a sigmoid truncated by some cut-off function), the exponential stability result could hold by making use of the Łojasiewicz inequality: for a compact set K ⊂ R d and two continuous, sub-analytic functions g, h :

K -→ R, if g -1 ({0}) ⊂ h -1 ({0}) then h(x) α c g(x)
for all x ∈ K, holds for some c > 0 and α ∈ N (see [START_REF] Lojasiewicz | Sur le problème de la division[END_REF] and also [Bierstone and Milman, 1988, Theorem 6.4]). This however remains an open problem. On the other hand, addressing analytically the (exponential) stability stipulated by the numerical experiments presented herein for non 2 -losses such as cross-entropy also remains open.

A pre-training algorithm. Note that the insight from Theorem 4.1 and the experiments performed for (4.7) could motivate the following algorithm, which trains a shallower network (i.e., a neural ODE with a shorter time horizon) at first, and then successively increases this time horizon in an iterative manner. One of the most distinguished characteristics of (4.2) is that the time-horizon T needed to get ε-close to any given target is in fact implicitly defined in the cost functional. At the level of ResNets, this means that the required number of layers needed to fit the data up to ε-error is given by the cost itself. The goal of the above algorithm is to take advantage of this artifact, and represents a greedy algorithm which uses only the number of layers strictly needed, thus avoiding unnecessary ones. We leave open the possible numerical analysis of the above algorithm (see [Goodfellow et al., 2016, Chapter 15] for some related literature on pre-training algorithms).

More general models. It would be of interest to directly investigate the appearance of the phenomena presented herein to neural ODEs with state of the art configurations, including convolutional, batch-normalization and max-pooling layers -we have solely focused our theoretical analysis to the basic settings. Further extensions could include the study of ResNet variants such as MomentumNets [START_REF] Sander | Momentum residual neural networks[END_REF], which are second order ODEs wherein ẋ is seen as a damping term, and mean field variants [START_REF] Weinan | A mean-field optimal control formulation of deep learning[END_REF].

Unsupervised learning. As discussed in the introduction, the neural ODE representation of deep supervised learning, due to the invertibility of the neural ODE flow, has seen fruitful applications in generative modeling via continuous normalizing flows (see [START_REF] Ruthotto | An introduction to deep generative modeling[END_REF]). In generative modeling and unsupervised learning, one aims to infer the probability distribution of the inputs { x i } i∈[N ] rather than, for instance, draw a decision boundary as in classification tasks via supervised learning. It would be of interest, in view of the existing applications, to investigate the potential use of the results presented in this work to the context of unsupervised learning. As c n -→ 0 along any subsequence as n → ∞ by virtue of the weak convergences of {w n } ∞ n=1 to w † and {b n } ∞ n=1 to b † , we deduce (A.1). Hence, Φ T sends bounded sequences in L 2 (0, T ; R du ) into strongly convergent (sub)sequences in C 0 ([0, T ]; R dx ) and is thus compact. Proof of (ii). Let {[w n , b n ]} ∞ n=1 ⊂ L 2 (0, T ; R du )∩BV(0, T ; R du ) be a bounded sequence in L 2 (0, T ; R du ) ∩ BV(0, T ; R du ). By the compactness of the embedding BV(0, T ; R du ) → L 1 (0, T ; R du ) (see [Ambrosio et al., 2000, Theorem 3.23]), there exists a pair w † , b † ∈ BV(0, T ; R du ) such that, along some subsequence, Arguing as above, we see that

x n (t)x † (t) c(σ) Using the fact that x n is bounded uniformly in n in C 0 ([0, T ]; R dx ) (this follows by applying a Grönwall argument) and the strong L 1 -convergences of the parameters, we conclude that c n -→ 0, whence (A.1) follows.

The proof of (iii) follows by arguing as in (i) and (ii), with an intermediate use of the Rellich-Kondrachov compactness theorem to ensure strong convergence of the parameters in L 2 -we omit the proof. This concludes the proof.

Proof of Lemma 7.3. We focus on proving (i) and split the proof in two steps. The proof of (ii) will follow by simply applying a Cauchy-Schwarz inequality at the conclusion of each of the two steps.

Step 1. Let us first suppose that t 1. By integrating (3.2) on the interval [0, t] ⊂ [0, 1] and using the fact that σ ∈ Lip(R) and σ(0) = 0, it may be seen that x(t)

x 0 + c(σ) x(t) w L 1 (0,1;R d×d ) + b L 1 (0,1;R d )

x 0 -x + c(σ)N 1 /2 max 1, C [w, b] L 1 (0,1;R du ) [w, b]
L 1 (0,1;R du ) .

Step 2. Now suppose that t ∈ (1, T ]. We first show that there exists a t * ∈ (t -1, t] such that x(t * )x xx L 2 (0,T ;R dx ) .

(A.3) This follows by a contradiction argument -indeed, suppose that x(t * )x > xx L 2 (0,T ;R dx ) . for all t * ∈ (t -1, t], then T 0 x(s)x 2 ds t t-1

x(s)x 2 ds > T 0

x(s)x 2 ds, which is a contradiction. Consequently, we know that there exists a t * ∈ (t -1, t] such that (A.3) holds. By integrating (3.2) in [t * , t] and using the fact that σ ∈ Lip(R) and σ(0) = 0, we see that

x(t) -x x(t * ) -x + c(σ) t t *
w(s) x(s)x + x w(s) + b(s) ds.

By virtue of Grönwall's inequality, it can then be seen that

N -1 /2 x(t) -x C x(t * ) -x + [w, b]
L 1 (0,T ;R du ) exp C w L 1 (0,T ;R d×d )

for some C = C (σ, x ) > 0. Using (A.3), we conclude the proof.

Proof of Lemma 7.4. We will make use of the notation u := [w, b] for simplicity.

Let us first suppose that T 1 + τ • . In this case, the arguments follow precisely those presented in the proof of Theorem 4.2. We just give a brief sketch of the main differences. Due to the controllability assumption, there exist u † ∈ L 2 (τ • , 1 + τ • ; R du )

  is either L 2 or H 1 . In Theorem 3.1, we show that when loss(•, •) and the affine map P : R d → R m are such that the minimum of E (equal to zero) is attained and when the activation function σ in (1.2) is 1-homogeneous, the training error E(x T (T )) of the collection x T = {x T,i } i∈[N ] of solutions to (1.2) corresponding to any solution [w T , b T ] to the minimization problem (1.3), is at most of the order O 1 T , whilst the optimal parameters [w T , b T ] converge, on a suitable time-scale, to a solution [w * , b * ] of inf [w,b]∈H(0,1;R du ) x i (•) solves (1.2) in [0

Figure 1 .

 1 Figure 1. Example 4.1: We depict a manifestation of the stability results of Theorem 4.1 for the state x T (t) (right) and the training error E(x T (t)) (left) over t ∈ [0, T ]. We observe that, after a finite time, the training error and trajectory remain at a steady configuration, so further times/layers could be discarded from training.

Figure 2 .

 2 Figure2. Example 4.1: A batch of training data (left) and the evolution of the corresponding trained trajectories x T,i (t) (right) in the phase plane. The learned flow is simple, with moderate variations, due to the exponentially small parameters. We refer to https://github.com/borjanG/dynamical.systems/blob/master/ videos/example4-1.mp4 for a movie of the evolution of the trajectories, where the stability phenomenon depicted in Figure1can be visualized.

Figure 3 .

 3 Figure 3. Example 4.1:The trained classifier on [-2, 2] 2 and evaluated on a batch of the test set. The simplicity of the learned flow ensures relatively satisfactory generalization as the shape of the dataset is learned adequately, and the test set is correctly classified.

Figure 4 .

 4 Figure 4. Example 4.2: The decay of the training error (left) and stabilization of the trained trajectories x(t) and {P x i (t)} i∈[N ] (right).

  Example 4.3 (Annuli in three colors). We consider a toy classification task with three labels, namely y i ∈ [3], with each label corresponding to a different color. The dataset consists of N = 3200 training samples and 800 test samples. We consider the crossentropy loss (3.10) in the training error in (4.7) and λ = 0.01, and we consider the neural ODE (3.3) with T = 15 (we use a forward Euler scheme to obtain a corresponding ResNet with fixed time-step equal to 1), where σ ≡ tanh. No augmentation of the initial8

Figure 5 .

 5 Figure 5. Example 4.2: A batch of training data (left) and the evolution of the corresponding trained trajectories x T,i (t) (right). See https://github.com/borjanG/dynamical.systems/blob/master/ videos/exampleA-1.mp4 for a movie of the evolution of the trajectories, where the stability phenomenon depicted in Figure 4 can be visualized.

Figure 6 .

 6 Figure 6. Example 4.2:The trained classifier on [-2, 2] 2 and its evaluation on a batch of the test data. The shape of the data is captured accurately.

Figure 7 .

 7 Figure 7. Example 4.3: The decay of the training error (left) and stabilization of the trained trajectories x(t) and {P x i (t)} i∈[N ] (right).

Figure 8 .

 8 Figure 8. Example 4.3: A batch of training data (left) and the evolution of the corresponding trained trajectories x T,i (t) (right). See https://github.com/borjanG/dynamical.systems/blob/master/ videos/example4-2.mp4 for a movie of the evolution of the trajectories, where the stability phenomenon depicted in Figure 7 can be visualized.

Figure 11 .

 11 Figure 11. Example 4.4: A batch of training data (left) and the evolution of the corresponding trained trajectories x T,i (t) (right). See https://github.com/borjanG/dynamical.systems/blob/master/ videos/example4-3.mp4 for a movie of the evolution of the trajectories, where the stability phenomenon depicted in Figure 10 can be visualized.

Figure 12 .

 12 Figure 12. Example 4.4:The trained classifier on [-2.5, 2.5] 2 and its evaluation on a batch of the test data. The shape of the dataset is captured accurately.

Figure 13 .

 13 Figure 13. Example 4.5: The decay of the training error (left) and stabilization of the trained trajectories x(t) and {P x i (t)} i∈[N ] (right).

Figure 10 .Figure 11 .

 1011 Figure 10. Example 4.4: The decay of the training error (left) and stabilization of the trained trajectories x(t) and {P x i (t)} i2[N ] (right).

Figure 14 .

 14 Figure 14. Example 4.5:We depict the evolution of two individual samples x i (t) ∈ R 784 mapped onto a 28 × 28 grid. We see that each trajectory stabilizes to some stationary configuration. The trained model tends to "diffuse" (in a caloric sense) the input signal ahead of classifying via the softmax applied to P x i (t) ∈ R 10 .

Figure 15 .Figure 16 .

 1516 Figure 15. Example 4.5: The validation error and accuracy over training epochs (experiments repeated 10 times); in this simplified dataset setting, generalization is not necessarily compromised due to the introduction of an integrated empirical risk.

Figure 18 .Figure 18 .

 1818 Figure18. Example A.2: The validation error and accuracy over training epochs (experiments repeated 10 times); as anticipated, generalization is not as good as for the simpler MNIST dataset. The lower accuracy with respect to state of the art configurations could also be due to the fact that we do not make use of convolutional layers.

  (m -1) exp -γ exp T α 2 .

  .22) By Lemma 7.5, the parameters [w T , b T ]| [τ•,T ] minimize the functional J τ•,T for System (7.20) with initial data x τ• = x T (τ • ), to which the solution is precisely x T | [τ•,T ] . Now, applying Lemma 7.4 (this time with with τ • as in (7.22), x τ• = x T (τ • ) and with r = x 0x , as C1 / √ τ < 1) in combination with (7.22), yields

  .23) which holds for all t ∈ [τ • , T ]. Since τ • < τ , (7.23) also holds for t ∈ [τ, T ].

Algorithm 1 :

 1 A pre-training algorithm. Result: T , [w T , b T ] fix ε > 0, T • > 0; j ← 1; Find [w jT• , b jT• ] and x jT• by solving inf [w,b]∈H k ((j-1)T•,jT•;R du ) = f [w(t), b(t)], x(t) in ((j -1)T • , jT • ) x((j -1)T • ) = x (j-1)T• . (8.2) [w T , b T ] | ((j-1)T•,jT•) ← [w jT• , b jT• ]; x j ← x jT• (jT • ); T ← jT • ; while E(x(jT • )) > ε do j ← j + 1; find [w jT• , b jT• ] and x jT• by solving (8.1) -(8.2); [w T , b T ] | ((j-1)T•,jT•) ← [w jT• , b jT• ]; x jT• ← x jT• (jT • ); T ← jT • end

  (τ ) w n (τ )w † (τ ) dτ + t 0 b n (τ )b † (τ ) dτ .Using Grönwall's inequality, Cauchy-Schwarz, and the boundedness of the L 2 -norm of {w n } ∞ n=1 by some constant M > 0 independent of t, we thence obtainx n (t)x † (t) c n exp c(σ)

[

  w n , b n ] ----→ n-→∞ w † , b † strongly in L 1 (0, T ; R du ). Of course, the same convergences thence hold for w n := diag N (w n ) to w † := diag N (w † ), as well as b n := [b n , . . . , b n ] to b † := [b † , . . . , b † ]. Let x † ∈ C 0 ([0, T ]; R dx ) be the unique solution to (3.2) associated to [w † , b † ] and the initial datum x 0 . Let us prove (A.1).

t 0 w

 0 n (τ )w † (τ ) x n (τ ) dτ + c(σ) t 0 w † (τ ) x n (τ )x † (τ ) dτ + c(σ) t 0 b n (τ )b † (τ ) dτ,where c(σ) > 0 is the Lipschitz constant of σ, and thus, using Grönwall's inequality, we obtainx n (t)x † (t) c(σ)c n exp c(σ) τ )w † (τ ) x n (τ ) dτ + t 0 b n (τ )b † (τ ) dτ.

  for t ∈ [0, 1], where c(σ) > 0 denotes the Lipschitz constant of σ. By the Grönwall inequality, we then havex(t) x 0 + c(σ)N 1 /2 b L 1 (0,1;R d ) exp c(σ)N 1 /2 w L 1 (0,1;R d×d ) :=C [w,b] L 1 (0,1;R du ) (A.2) for t ∈ [0, 1]. Using (A.2), it may be seen thatx(t)x x 0x + c(σ)N 1 /2

  Figure 10. Example 4.4: The decay of the training error (left) and stabilization of the trained trajectories x(t) and {P x

				Decay of training error						Stability of norms
	1.0							E(x(t))	160 180
	0.8								140
									120
	0.6								100
	0.4								80
									60
	0.2								40
	0.0								20							|x(t)| 2 |P x(t)| 2
	0	2	4	6	8	10	12	14	0	2	4	6	8	10	12	14
					t (layers)								t (layers)
	train	Generalization outside training data		2.00
	test								1.75
									1.50
									1.25		Figure 9. Example 4.3:
	x 2							1.00		The trained classifier
									0.75		on [-2.5, 2.5] 2 and its
									0.50		evaluation on a batch
					x 1				0.25 0.00		of the test data. The shape of the dataset is captured accurately.

T,i (t) (right). See https://github.com/borjanG/dynamical.systems/blob/master/ videos/example4-2.mp4 for a movie of the evolution of the trajectories, where the stability phenomenon depicted in Figure

7

can be visualized. grayscale, 28 × 28 image of a handwritten digit, and thus an element of R 784 ; the dataset has 10 labels: y i ∈ [10]. We consider a similar setup as in Example 4.4 -the model is i (t)} i∈[N ] (right).

  the training error in (4.7) and only L 2 -regularization of the parameters, with T = 20. We emphasize that we do not use any convolutional layers nor other commonly used operations in image classification (e.g. batch normalization, max-pooling) in the underlying ResNet architecture, and we solely concentrate on basic matrix weights. The output layer is parametrized by P x = p 1 x + p 2 , where p 1 ∈ R 10×784 , p 2 ∈ R 10 are part of the trainable parameters. We show the results of the experiments in Figure13-Figure14.

	0.12			Decay of training error		E(x(t))	450			Stability of norms
	0.10
										400
	0.08									350
	0.06									300
	0.04									250
	0.00 0.02									200 150							|x(t)| 2 |P x(t)| 2
	0.0	2.5	5.0	7.5	10.0	12.5	15.0	17.5	20.0	0.0	2.5	5.0	7.5	10.0	12.5	15.0	17.5	20.0
					t (layers)									t (layers)
										train		Generalization outside training data		1.05
										test								0.90
																		0.75
																		0.60
										2
										x
																		0.45
																		0.30
																		0.15
														x 1				0.00

parametrized as (2.4) -(2.7) (we use a forward Euler scheme to obtain a corresponding ResNet with fixed time-step) where d hid = 32 and σ ≡ tanh, we consider cross-entropy loss in

  and thus x k ∈ R d k for k ∈ {0, . . . , N layers }, where {d k }

	N layers
	k=0

Here H 1 (0, T ; R du ) denotes the Sobolev space of square integrable functions from (0, T ) to R du with square integrable weak derivatives (see Section 1.4). We consider H 1 -regularization in the setting of (1.2) to ensure the existence of minimizers, see Remark 1.

CARLOS ESTEVE-YAGÜE, BORJAN GESHKOVSKI, DARIO PIGHIN, AND ENRIQUE ZUAZUA

Acknowledgments. B.G. acknowledges Daniel Tenbrinck and Lukas Pflug (FAU Erlangen-Nürnberg) for discussions on the foundations of neural networks and non-local equations. Funding: B.G. and E.Z. have received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No.765579-ConFlex. D.P., C.E. and E.Z. have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement NO. 694126-DyCon). The work of E. Z. has been supported by the Alexander von Humboldt-Professorship program, the Transregio 154 Project "Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks" of the German DFG, grant MTM2017-92996-C2-1-R COSNET of MINECO (Spain) and by the Air Force Office of Scientific Research (AFOSR) under Award NO. FA9550-18-1-0242.

Combining (7.7) and (7.6), we deduce that

L 2 (0,1;R du ) (7.8) for any T > 0. Now by integrating (3.3), and using the fact that σ is globally Lipschitz continuous with constant C(σ) > 0 and satisfies σ(0) = 0, for any t ∈ [0, T ] we have

By using the Grönwall inequality, we obtain

whereas by Cauchy-Schwarz, it follows that

At this point, employing (7.8), we deduce

Since u 1 is independent of T , we conclude that the set {x T (T )} T >0 is bounded. Whence, there exists a sequence {T n } ∞ n=1 with T n > 0 and T n -→ ∞ as n -→ ∞ and some

(7.9)

Since E (x Tn (T n )) -→ 0 as n -→ ∞ by (7.4), by continuity of E, we have E(x • ) = 0. This concludes the proof of (3.6).

Part 3. We now address the third statement of the theorem. To this end, we will first show that the sequence {u n } ∞ n=1 defined in the statement is bounded in L 2 (0, 1; R du ). Let u † ∈ L 2 (0, 1; R du ) be any solution to

We prove this claim by contradiction. Indeed, assume that we had

for some n 1.

We consider

(1) by Lemma 7.1. On another hand, by assumption we have E(x † (1)) = 0. It then follows Appendix A. Auxiliary proofs Proof of Lemma 7.2. The proof is mostly identical in for all of the three items, the difference being an underlying compact embedding.

By the Banach-Alaoglu theorem, there exists a pair w † , b † ∈ L 2 (0, T ; R du ) such that, along some subsequence,

Of course, the same convergences thence hold for

] and the initial datum x 0 . Let us prove that

along the aforementioned subsequence. Take an arbitrary t ∈ [0, T ]. Note that

Hence, using the fact that σ is globally Lipschitz with constant c(σ) > 0,

By a Grönwall argument we see that for t ∈

and we denote by x aux (•) the associated solution to (7.20). We note that x aux (t) = x for t ∈ [1 + τ • , T ] and thus E(x aux (T )) = 0. Just as for the proof of Theorem 4.2, by using J τ•,T (u T ) J τ•,T (u aux ), the fact that E(x aux (T )) = 0, E 0, (A.4), (4.4), and recalling Lemma 7.3, we may conclude the proof when

we may combine this inequality with Lemma 7.3 to conclude the proof.

Proof of Lemma 7.5. We argue by contradiction. Suppose that there exist parameters w

Denoting by x aux the associated neural ODE trajectory on [0, T ] with x aux (0) = x 0 , we see that x aux (T ) = x † (T ) and so E(x aux (T )) = E(x † (T )). Consequently, . We argue by contradiction. Suppose that f (t) X > η(τ ) for all t ∈ [a, a + τ ).

which is a contradiction.