Joubine Aghili

Olga Mula

An optimal control framework for adaptive neural ODEs

An optimal control framework for adaptive neural ODEs

Joubine Aghili, Olga Mula

Introduction

Context

Neural networks produce structured parametric families of functions that have been studied for at least 70 years (see [START_REF] Hebb | The Organization of Behavior: A Neuropsychological Theory[END_REF][START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF]). It is however only in the last decade that their popularity has surged. Thanks to the increase of computing power, and the development of easy-to-use computing tools for optimization and automatic differentiation, statistical learning of neural networks has produced state-of-the-art performance in a large variety of machine learning problems, from computer vision [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] (e.g. self-driving cars, X-rays diagnosis,...) to natural language processing [START_REF] Wu | Google's neural machine translation system: Bridging the gap between human and machine translation[END_REF] (e.g. Google translate, DeepL Translator,...) and reinforcement learning (e.g. superhuman performance at Go [START_REF] Silver | Mastering the game of go with deep neural networks and tree search[END_REF]). Despite this great empirical success, neural networks are not entirely well-understood and there is a pressing need to provide:

• solid mathematical foundations to understand their approximation power and why they may outperform other classes of functions, • algorithms to find systematically an appropriate architecture for each problem and to train the network in order to deliver in practice the approximation capabilities predicted by the theory, and to guarantee robustness over generalization errors. On the first point, significant advances have been made recently on the approximation power of neural networks (see, e.g., [START_REF] Yarotsky | Error bounds for approximations with deep ReLU networks[END_REF][START_REF] Daubechies | Nonlinear approximation and (deep) relu networks[END_REF][START_REF] Grohs | Deep neural network approximation theory[END_REF][START_REF] Gühring | Error bounds for approximations with deep ReLU neural networks in W s,p norms[END_REF] for a selection of rigorous results). In particular, a few recent works have given theoretical evidence on the advantages of using deep versus shallow architectures for the approximation of certain relevant families of functions (see, e.g., [START_REF] Daubechies | Nonlinear approximation and (deep) relu networks[END_REF][START_REF] Telgarsky | Representation benefits of deep feedforward networks[END_REF]). However, so far the obtained results do not seem to be informative on how to address the second point above, that is, how to build algorithms that discover automatically the right architecture for each problem, and that may allow to benefit from the theoretically high approximation properties of (deep) neural networks. Even though there are many contributions using evolutionary algorithms, e.g. [START_REF] Angeline | An evolutionary algorithm that constructs recurrent neural networks[END_REF][START_REF] Sun | Automatically designing CNN architectures using the genetic algorithm for image classification[END_REF], the theoretical foundations of adaptive neural networks remains unclear. As a result, their training remains a key open issue subject to very active research. The present work is a contribution in this direction from the perspective of neural ODEs [START_REF] Chen | Neural ordinary differential equations[END_REF].

The most commonly applied training methods are based on the stochastic gradient descent ( [START_REF] Robbins | A stochastic approximation method[END_REF][START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]) and its variants (see, e.g., [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF][START_REF] Zeiler | Adadelta: an adaptive learning rate method[END_REF]). It has the advantage of being easily implementable but it is difficult to tune properly in practice, especially in problems involving large data sets and using deep neural networks with many layers and coefficients which are very prone to over-parametrization issues. The difficulty in training deep networks raises the question regarding the benefits of applying a shallow-to-deep adaptive strategy in the network architecture in order perform the training in a more robust and rapid manner. The underlying intuition is that shallow networks give poor approximation properties but converge faster than deeper networks which, on the contrary, have higher approximation power.

There is a large body of emerging works which explore numerically the potential of training deep neural networks with shallow-to-deep strategies (see, e.g., [START_REF] Chen | Net2net: Accelerating learning via knowledge transfer[END_REF][START_REF] Wei | Network morphism[END_REF][START_REF] Li | Learning without forgetting[END_REF][START_REF] Cortes | Adanet: Adaptive structural learning of artificial neural networks[END_REF][START_REF] Elsken | Neural architecture search: A survey[END_REF][START_REF] Yang | Adaptive deep models for incremental learning: Considering capacity scalability and sustainability[END_REF]). In this work, we present an abstract setting and convergence results which should be understood as one possible theoretical justification for this type of adaptive approach. In order to understand the potential gain, it is necessary to define a notion of the continuous underlying objects that are approximated. This is the reason why we have chosen to work from a perspective which connects the task of learning with neural ODEs with optimal control and dynamical ODE systems. In this view, there is a notion of a continuous neural network which is described by a continuous time-dependent ODE system. Its discretization fixes its architecture. For instance, ResNet (see [START_REF] He | Deep residual learning for image recognition[END_REF]) can be regarded as an Explicit Euler scheme for the solution of certain dynamics. The weights involved in the ODE are seen as controls and correspond to the hyperparameters of the network which have to be optimized. The process of learning these parameters can then be recast as an optimal control problem over the admissible controls where the cost function is the empirical risk (sometimes with an extra regularization term). Necessary optimality conditions are then easily formulated through the Pontryagin Maximum Principle (PMP, [START_REF] Boltyanskii | The theory of optimal processes. i. the maximum principle[END_REF][START_REF] Pontryagin | Mathematical Theory of Optimal Processes[END_REF]).

Contribution of the present work

The optimal control point of view is appealing since it gives access to a fully continuous description of neural networks and of the underlying statistical learning process. The practical implementation involves two discretization errors:

1. Sampling error: The cost function of the fully continuous optimal control problem is the average of a loss function, sometimes also called risk. However, in practice, the average loss is replaced by an empirical mean which is built from N samples. So it is necessary to assess by how much the minimum of the sampled mean deviates from the continuous underlying problem. We analyze this point in Theorem 7.2. 2. Time-integration error: To solve the optimal control problem on the empirical mean of the loss, one solves the resulting PMP with Newtontype iterations involving the solution of forward and backward dynamical systems with certain time integration schemes (which fix the neural network architecture). The usual approach is to first fix a discretization and then solve the discrete version of the PMP. However, the discretization errors accumulate at each iteration and make us deviate from the timecontinuous version. It is thus necessary to examine how much the final output deviates from the time-continuous one. In this context, the main contributions of the present work are:

1. Derivation of a full error analysis:

• We make a theoretical analysis on the interplay between the above two errors. The impact of the sampling error is analyzed in Theorem 7.2 and the discretization error in Theorem 6.2. From these results, it follows that to approximate the minimum of the underlying fully continuous problem at a certain accuracy, we not only need a minimal number of training samples, but we also need to solve the control problem on the empirical risk with some minimal accuracy. This requires, in turn, to use time-adaptive techniques. In our analysis, both errors are additive and can actually be optimally balanced. The final result is given in Corolary 7.3 and it is formulated in probability. • To solve the control problem on the empirical risk at a given target accuracy, we show that we have to compute the forward and backward propagations at each step of the Newton-type algorithms at increasing accuracy. As our theoretical proof of Theorem 6.2 will illustrate, tightening the accuracy of the time integration is necessary to guarantee convergence to the exact continuous problem and justifies the coarse-to-fine, resp. shallow-to-deep, strategy. This yields neural networks with depths that are automatically adapted. 2. Practical coarse-to-fine adaptive algorithm:

• We next use the theory to derive actionable criteria to build timeintegration schemes that are adaptively refined across the iterations. This has the advantage of reducing computational cost since early iterations use looser tolerances, thus avoiding unnecessary work due to oversolving, while later iterations use tighter tolerances to deliver accuracy. The coarse-to-fine strategy could also be understood as a dynamic regularization which helps to prevent over-parametrization. • We illustrate the behavior of the algorithm in several benchmark examples.

Potential impact, limitations and extensions

We believe the results of this article can contribute to the following key topics in deep learning:

Related works

To the best of our knowledge, the connection between deep learning, dynamical systems and optimal control is not new. It can be traced back at least to the works of LeCun and Pineda in the 1980's, in which the idea of back-propagation is connected to the adjoint variable arising in optimal control (see [START_REF] Lecun | A theoretical framework for back-propagation[END_REF][START_REF] Pineda | Generalization of back propagation to recurrent and higher order neural networks[END_REF]).

The approach has gained a lot of attention in recent years following works on neural ODEs such as [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Weinan | A proposal on machine learning via dynamical systems[END_REF][START_REF] Haber | Stable architectures for deep neural networks[END_REF]. We refer to these works, and also to [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][START_REF] Vialard | A Shooting Formulation of Deep Learning[END_REF][START_REF] Ayyubi | Progressive growing of neural odes[END_REF] for some selected references.

In the field of neural ODEs, there are by now many experimental studies exploring the effects of different time-stepping techniques (see e.g. [START_REF] Lu | Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations[END_REF][START_REF] He | Ode-inspired network design for single image super-resolution[END_REF]), and their effects on stability and robustness (see, e.g., [START_REF] Haber | Stable architectures for deep neural networks[END_REF][START_REF] Chang | Reversible architectures for arbitrarily deep residual neural networks[END_REF]). The topic of adaptivity has also been explored. In works such as [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Zhang | Anodev2: A coupled neural ode framework[END_REF][START_REF] Zhuang | Adaptive checkpoint adjoint method for gradient estimation in neural ode[END_REF], high order adaptive methods which automatically change the time step by estimating local truncation error have been developed . We also refer to [START_REF] Benning | Deep learning as optimal control problems: Models and numerical methods[END_REF] where the time step sizes are included as an additional parameter to be optimized. In all these works the number of layers is fixed during the training iterations. This idea is part of our work but we also understand adaptivity as the idea of increasing progressively the number of steps. In fact, given the present state of the art, our contribution should be undertsood as an abstract, theoretical framework that rigorously justifies convergence to the continuous problem, rather than a very concrete numerical algorithm. Our numerical implementation does not provide state of the art results but it is intended to reproduce as faithfully as possible the theoretical setting.

Outline of the paper

In section 2, we define precisely what we understand by deep learning and we recall its connection with optimal control. Section 3 gives the optimal control setting which we use in our subsequent developments. In particular, we formulate the fully continuous optimal control problem over the expectation of the loss and its sampled version involving an empirical mean of the loss. Section 4 recalls the Pontryagin's Maximum Principle for the sampled problem. We also recall in this section the Extended Method of Successive Approximations (E-MSA), originally introduced in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF], to solve the sampled problem. The algorithm is iterative, and it requires to compute at each iteration a forward and a backward time propagation followed by a certain maximization over the control variables. In its original formulation, it is assumed that these steps are performed exactly. In Section 5, we formulate our main algorithm, which is an inexact version of the E-MSA that we call Adaptive MSA (A-MSA). It is based on approximately realizing each propagation and maximization within a certain accuracy. To ensure convergence to the exact solution, the accuracy needs to be tightened across the iterations. In section 6, we prove that A-MSA converges to a time-continuous solution of the sampled PMP. In section 7, we connect the solution of the sampled problem with the continuous one involving the exact average of the loss function. For this, we use results obtained in [START_REF] Han | A mean-field optimal control formulation of deep learning[END_REF]. This last step allows to connect how much the solution of A-MSA deviates from the fully continuous one in terms of the number of samples and the accuracy of the discretization. In section 8, we give guidelines to implement A-MSA in practice, and section 9 illustrates the performance of the algorithm in numerical examples. We conclude the paper in Section 10.

From deep learning to optimal control

Statistical Learning Problems

We consider the following statistical learning problem: Assume that we are given a domain set X ⊆ R n and a label set Y ⊆ R k , with n, k ∈ N. Further assume that there exists an unknown probability distribution µ on X × Y .

Given a loss function L : Y × Y → R + , the goal of the statistical learning problem is to find a function v : X → Y , which we will call prediction rule, from a hypothesis class V ⊂ {v : X → Y } such that the expected loss

J (v) := E (x,y)∼µ L(v(x), y)
is minimized over V. In other words, the task is to find

v * ∈ arg min v∈V J (v).
In general, the probability distribution µ is unknown and we are only given a set S N of N ∈ N training samples

S N := {(x i , y i )} N i=1 .
The most common choice is then to consider the uniform distribution µ N :

X × Y → R + , µ N (x, y) = 1 N N i=1 δ (xi,yi) (x, y)
which yields the so-called empirical loss

J N (v) := E (x,y)∼µ N L(v(x), y) = 1 N N i=1 L(v(x i ), y i )
and one finds v N , an approximation of v, by minimizing over it,

v N ∈ arg min v∈V J N (v).
(2.1)

In the following, the above optimization step will be called the learning procedure.

Neural Network Architectures and Deep Learning

While there exists numerous different architectures, we could simplify the discussion and say that neural networks are a class of functions with the basic general form

v : R n → R k , x → W L σ(W L-1 σ(. . . (σ(W 1 (x)))), (2.2)
where:

• L ∈ N is the depth, i.e., the number of layers of the neural network,

• W ℓ : R N ℓ-1 → R N ℓ are affine maps for ℓ = 1, . . . , L. For consistency, we must set N 0 = n and N L = k but the rest of the dimensions N ℓ ∈ N can be freely chosen. We have for all

x ∈ R N ℓ-1 , W ℓ (x) = A ℓ (x) + b ℓ for a matrix A ℓ ∈ R N ℓ ×N ℓ-1 and a vector b ℓ ∈ R N ℓ .
• σ : R → R is a nonlinear activation function which is applied coordinatewise in (2.2). Some popular choices are the ReLU function, σ(x) = max(0, x), or the hyperbolic tangent, σ(x) = tanh(x). Deep learning describes the range of learning procedures to solve statistical learning problems where the hypothesis class V is taken to be a set of neural networks. One usually works with the class of neural networks with given depth L, activation function σ and fixed dimensions N 1 , . . . , N L-1 for the affine mappings,

N N (L, σ, N 1 , . . . , N L-1 ) := {v : X → Y : v(x) = W L σ(W L-1 σ(. . . (σ(W 1 (x))))}.
For this class, the task is to solve the empirical risk problem (2.1) associated to it.

Continuous Formulation as an Optimal Control Problem

Functions of the type (2.2) can be built by repeated composition of parametrizable functions

ϕ ℓ : R N ℓ-1 → R N ℓ , x → ϕ ℓ (x) = σ(W ℓ (x)), i = 1, . . . , L -1,
followed by an affine step, that is,

v = W L • ϕ L-1 • • • • • ϕ 1 , ∀v ∈ N N (L, σ, N 1 , . . . , N L-1 ).
This simple observation yields to numerous other possible architectures. The most relevant instance triggering the connection with optimal control are the so-called Residual Neural Networks (ResNet, [START_REF] He | Deep residual learning for image recognition[END_REF]). They correspond to the choice

ϕ ℓ : R N ℓ-1 → R N ℓ , x → ϕ ℓ (x) = x + hσ(W ℓ (x)), i = 1, . . . , L -1,
for a certain parameter1 h > 0.

If we now set

N L-1 = • • • = N 1 = N 0 = n,
the application of ϕ ℓ can be interpreted as an Explicit Euler step from time t ℓ-1 = (ℓ -1)h to time

t ℓ = t ℓ-1 + h of the dynamics ẋt = f (x t , θ t ), with x 0 = x, where f (x t , θ t ) = σ(W t (x t )), with W t (x) = A t x + b t ,
Network Architecture Associated ODE discretization ResNet [START_REF] He | Deep residual learning for image recognition[END_REF], RevNet [START_REF] Gomez | The reversible residual network: Backpropagation without storing activations[END_REF] Forward Euler Polynet [START_REF] Zhang | Polynet: A pursuit of structural diversity in very deep networks[END_REF] Approximation of Backward Euler FractalNet [START_REF] Larsson | Fractalnet: Ultra-deep neural networks without residuals[END_REF] Runge-Kutta order 2

Table 1: Some neural network classes and their associated ODE discretization. and θ t gathers the parameters upon which the dynamics depend. In our case,

θ t = (A t , b t ).
Note that the matrices A and b are now time-dependent, hence the t-subscript.

In this view, ResNet functions can be interpreted as the output of performing L -1 time steps of size h of the above dynamics, followed by an affine operation W L : R n → R k . This last operation is important since it maps the final output from the domain set in R n to the label set in R k . In fact, it can be replaced by any linear or nonlinear mapping g : R n → R k . Note also that h can be included as a parameter of the model class and that one can optimize over it as in [START_REF] Benning | Deep learning as optimal control problems: Models and numerical methods[END_REF].

Following similar lines, one can view certain classes of neural networks as discretizations of an underlying continuous dynamical system. Table 1 lists some popular classes and the associated numerical scheme (see, e.g., [START_REF] Lu | Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations[END_REF] for more details on these connections). We emphasize however that the connection with optimal control comes at the price of imposing that the input and intralayer dimensions are equal (N L-1 = • • • = N 1 = N 0 ), thus it cannot be made for all classes of neural networks.

Optimal control setting

We next define the notation for the optimal control setting which we will use in the rest of the paper. The notation is kept as consistent as possible with section 2 in order to further enhance the similarity between the usual deep learning approach and the present one.

Like in section 2, we assume that we are given a domain set X ⊂ R n and a label set Y ⊂ R k , with n, k ∈ N, and that there exists an unknown probability distribution µ on X × Y representing the distribution of the input-target pairs (x, y). Consider now a set of admissible controls or training weights Θ ⊆ R m . Usually, we set Θ = R m in deep learning but the present methodology allows to consider constraints we will denote Θ ⊆ [U min , U max ] m . Fix T > 0 and let f (feed-forward dynamics), Φ (terminal loss function) and R (regularizer) be functions

f : R n × Θ → R n , Φ : R n × R k → R, R : Θ → R.
Let L ∞ ([0, T ], Θ) be the set of essentially bounded measurable controls. In the following, we use bold-faced letters for path-time quantities. For example, θ := {θ t : 0 ≤ t ≤ T } for any θ ∈ L ∞ ([0, T ], Θ). For every control θ ∈ L ∞ ([0, T ], Θ) and every value of the random variable x ∈ X, we define the state dynamics u θ,x := {u θ,x t : 0 ≤ t ≤ T } as the solution to the ordinary differential equation (ODE),

uθ,x t = f (u θ,x t , θ t ), ∀t ∈ (0, T ) u θ,x 0 = x (3.1)
The ODE is stochastic and its only source of randomness is the initial condition x.

With this notation, the deep learning optimization problem can be posed as the optimal control problem of finding

J * = inf θ∈L ∞ ([0,T ],Θ) J (θ), subject to (3.1), (3.2)
where

J (θ) := E (x,y)∼µ [Loss(x, y, θ)] ,
and for any input-target pair (x, y) ∈ R n × R k , the loss function is defined as

Loss(x, y, θ) := Φ(u θ,x T , y) + T 0 R(θ t )dt
Note that the regularizer R could in general also depend on the state u and Φ is the actual loss function upon which we want to act. It plays the same role as the loss function L of section 2.1. We can relate them by introducing the mapping g : R n → R k from section 2.3, and setting

Φ(u θ,x T , y) = L(g(u θ,x T ), y).
Like in the setting of section 2.1, we are only given a set of N samples S N = {(x i , y i )} N i=1 and the probability distribution µ is not exactly known. For a given control θ ∈ L ∞ ([0, T ], Θ), each sample follows the dynamics

uθ,i t = f (u θ,i t , θ t ), ∀t ∈ (0, T ) u θ,i 0 = x i (3.3)
for i = 1, . . . , N . We then perform empirical risk minimization taking, e.g., the uniform distribution

µ N = 1 N N i=1 δ (xi,yi)
. The optimal control problem becomes

J * S N = inf θ∈L ∞ ([0,T ],Θ) J S N (θ), subject to (3.3), (3.4)
where

J S N (θ) := E (x,y)∼µ N Loss(x, y, θ) = 1 N N i=1 Φ(u θ,i T , y i ) + T 0 R(θ t )dt
Note that the solutions of the sampled optimal control problem (3.4) depend on the sample set S N and are therefore random variables. However, for S N fixed, problem (3.4) is deterministic and can thus be solved with deterministic optimal control techniques. In what follows, we adopt this viewpoint first to solve problem (3.4) at any target accuracy. By this we mean the following: for any given target accuracy ε > 0, we prove in Theorem 6.2 that it is possible to numerically compute a control θ ε such that

J S N (θ ε ) -J * S N ≤ ε.
We then come back to the probabilistic point of view and prove in Corollary 7.3 that, provided that the number N of samples is sufficiently large, there exists a constant C > 0 such that

J S N (θ ε ) -J * ≤ Cε.
with high probabillity. The derivation of the latter bound relies on the results of [START_REF] Han | A mean-field optimal control formulation of deep learning[END_REF].

Pontryagin's Maximum Principle and Method of Successive Approximations

In this section, we focus on the sampled optimal control problem (3.4). We first recall the necessary optimality conditions, usually known as the Pontryangin's Maximum Principle. We next introduce the main algorithm to solve the PMP which will be the starting point for our subsequent developments. In the following, the Euclidean norm of any vector a ∈ R n is denoted by ∥a∥ and the scalar product with any other vector b ∈ R n is a • b.

Pontryangin's Maximum Principle

We define the Hamiltoninan H :

[0, T ] × R n × R n × Θ → R as H(t, u, p, θ) := p • f (u, θ) -R(θ). (4.1)
The following classical Pontryagin's Maximum Principle (PMP) gives necessary optimality conditions to problem (3.4).

Theorem 4.1. Let θ * ∈ L ∞ ([0, T ], Θ) be an optimal control to (3.4). For i = 1, . . . , N , let u θ * ,i be the associated state process with initial condition x i .

There exists a co-state process

p θ * ,i ∈ L ∞ ([0, T ], R n ) such that uθ * ,i t = f (u θ * ,i t , θ * t ), u θ * ,i 0 = x i , ṗθ * ,i t = -∇ u H(t, u θ * ,i t , p θ * ,i t , θ * t ), p θ * ,i T = -∇ u Φ(u θ * ,i T , y i ),
and, for each t ∈ [0, T ],

1

N N i=1 H(t, u θ * ,i t , p θ * ,i t , θ * t ) ≥ 1 N N i=1 H(t, u θ * ,i t , p θ * ,i t , θ), ∀θ ∈ Θ. (4.2)
The proof of this theorem and its variants can be found in any optimal control theory reference (see, e.g., [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF][START_REF] Clarke | The maximum principle in optimal control, then and now[END_REF][START_REF] Athans | Optimal Control: an Introduction to the Theory and Its Applications[END_REF]). We omitted the case of an abnormal multiplier, which will not be considered in the following.

We emphasize that the PMP is only a necessary condition, so there can be cases where the solutions to the PMP are not global optima for (3.4). Nevertheless, in practice the PMP often gives good solution candidates, and when certain convexity assumptions are satisfied the PMP becomes sufficient (see [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF]). In the next section, we discuss the numerical methods that we take as a starting point to solve the PMP.

MSA and Extended-MSA

A classical algorithm to find θ * and the corresponding forward and co-state dynamics u θ * ,i and p θ * ,i is the Method of Successive Approximations (see [START_REF] Chernousko | Method of successive approximations for solution of optimal control problems[END_REF]). It is a fixed-point method based on the following steps. Starting from an initial guess of the optimal control θ 0 , for each k ≥ 0 we first solve for i = 1, . . . , N the forward dynamics

uθ k ,i t = ∇ p H(t, u θ k ,i t , p θ k ,i t , θ k t ) = f (u θ k ,i t , θ k t ), u θ k ,i 0 = x i . (4.3)
The dynamics u θ k ,i allows us to compute the backward dynamics

ṗθ k ,i t = -∇ u H(t, u θ k ,i t , p θ k ,i t , θ k t ), p θ k ,i T = -∇ u Φ(u θ k ,i T , y i ). (4.4)
Finally, we update the control by using the maximization condition (4.2),

θ k+1 t ∈ arg max θ∈Θ 1 N N i=1 H(t, u θ k ,i t , p θ k ,i t , θ), ∀t ∈ [0, T ]. (4.5)
In this form, MSA converges only locally when it is initialized with a starting control guess that is sufficiently close to an optimal control θ * . To overcome this limitation, an Extended MSA algorithm (E-MSA) based on an augmented Lagrangian strategy has been introduced in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF]. This algorithm is the starting point for our subsequent developments. It works as follows: fix some ρ > 0 and define the augmented Halmitonian

H(t, u, p, θ, v, q) := H(t, u, p, θ) - ρ 2 ∥v -f (u, θ)∥ 2 - ρ 2 ∥q + ∇ u H(t, u, p, θ)∥ 2 . (4.6)
We can now formulate an extended PMP based on H. Proposition 4.2 (Extended PMP, see [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF]). Let θ * ∈ L ∞ ([0, T ], Θ) be an optimal control to (3.4). For i = 1, . . . , N , let u θ * ,i be the associated state process with initial condition x i . There exists a co-state process

p θ * ,i ∈ L ∞ ([0, T ], R n ) such that uθ * ,i t = ∇ p H(t, u θ * ,i t , p θ * ,i t , θ * t ), u θ * ,i 0 = x i , ṗθ * ,i t = -∇ u H(t, u θ * ,i t , p θ * ,i t , θ * t ), p θ * ,i T = -∇ u Φ(u θ * ,i T , y i ),
and, for each t ∈ [0, T ],

1

N N i=1 H(t, u θ * ,i t , p θ * ,i t , θ * t ) ≥ 1 N N i=1 H(t, u θ * ,i t , p θ * ,i t , θ), ∀θ ∈ Θ.
The E-MSA algorithm consists in applying the MSA algorithm with the augmented Hamiltonian H instead of with H. Since ∇ p H = ∇ p H = f and ∇ u H = ∇ u H, steps (4.3) and (4.4) remain the same as in MSA and the maximation step (4.5) is replaced by

θ k+1 t ∈ arg max θ∈Θ 1 N N i=1 H(t, u θ k ,i t , p θ k ,i t , θ), ∀t ∈ [0, T ]. (4.7)
It has been proven in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF] that if the parameter ρ is taken sufficiently large, the E-MSA algorithm converges to the set of solution of the extended PMP for any initial guess of the control θ 0 . However, this scheme cannot be realized in practice without discretizing. This introduces a perturbation of the timecontinuous formulation of the algorithm since the resulting trajectories will be approximations of the continuous one. To guarantee convergence to the exact continuous solution, it is necessary to identify suitable approximation error tolerances that still guarantee convergence to the exact solution. This motivates to introduce an Adaptive MSA scheme that we describe in the next section. We prove that if the forward and backward propagations are solved with increasing accuracy at each step, then the algorithm converges towards a continuous solution of the Extended PMP.

A-MSA: Adaptive Method of Successive Approximations

The algorithm requires defining forward and backward time-integration schemes, which we introduce in section 5.1. We next present the algorithm and how predictions are made in sections 5.2 and 5.3. We end up by proving convergence in section 6.

Routines solve_fwd and solve_bwd

In the following, we work with time-integration routines solve_fwd and solve_bwd which can be cast in the abstract framework that we next describe. Relevant particular instances of it will be continuous Galerkin (cG) or Runge-Kutta collocation (RK-C) methods.

Setting:

We describe the framework in the case of the forward dynamics

ut = f (u t ), ∀t ∈ (0, T ) u 0 = x (5.1)
for a generic Lipschitz continuous function f : R n → R n . The backward dynamics can be deduced similarly. Let 0 = t 0 < t 1 < • • • < t L = T be a partition of [0, T ] and let T ℓ := (t ℓ-1 , t ℓ ], and h ℓ := t ℓ -t ℓ-1 . For q ∈ N, let V T q be the space of continuous functions that are piecewise polynomials of degree q in the time mesh T = ∪ L ℓ=1 T ℓ , i.e.,

V T q := {v ∈ C 0 ([0, T ]; R n ) : v| T ℓ ∈ P (q) (T ℓ ), ℓ = 1, . . . , L},
where

P (q) (T ℓ ) := {v ∈ C 0 (T ℓ ; R n ) : ∀ t ∈ T ℓ , v(t) = q j=0 t j v j , v j ∈ R n }
is the space of polynomials of degree q in the interval T ℓ . Introducing a projection operator

Π (q) : C 0 ([0, T ]; R n ) → V T q
the time discrete approximation U to the solution u is defined as follows: we seek U ∈ V T q satisfying the initial condition U (0) = x as well as

U ′ t = Π (q-1) f (U t ), ∀t ∈ T ℓ , (5.2)
for ℓ = 1, . . . , L. Note that since both terms belong to P (q-1) (T ℓ ), (5.2) admits a Galerkin formulation

T ℓ v • U ′ t dt = T ℓ v • Π (q-1) f (U t ) dt, ∀v ∈ P (q) (T ℓ ). (5.3)
In the following, we use mainly (5.2) but (5.3) is of interest since it connects (RK-C) methods with (cG) methods. It is proven in [START_REF] Akrivis | Optimal order a posteriori error estimates for a class of runge-kutta and galerkin methods[END_REF], that the continuous Galerkin method corresponds to the choice Π (q-1) := P (q-1) , with P (q-1) denoting the (local) L 2 orthogonal projection onto P (q-1) (T ℓ ) for each ℓ. Runge-Kutta collocation methods with pairwise distinct nodes in T ℓ can be obtained by choosing Π (q-1) := I (q-1) with I (q-1) denoting the interpolation operator by elements of P (q) (T ℓ ) at the nodes

t ℓ-1 + α i (t ℓ -t ℓ-1 ), i = 1, . . . , q, ℓ = 1, . . . , L, with appropriate weights 0 ≤ α 1 < • • • < α q ≤ 1.
All RK-C methods with pairwise distinct nodes in [0, 1] can be obtained by applying appropriate numerical quadrature to continuous Galerkin methods.

L 2 error estimation: One simple way of estimating the error between u and U is via the following Grönwall inequality. Denoting e t = u t -U t , we have

ėt = f (u t ) -Π (q-1) f (U t ) = (f (u t ) -f (U t )) + f (U t ) -Π (q-1) f (U t )
Multiplying by e t , we get

1 2 d dt e 2 t = e t (f (u t ) -f (U t )) + e t f (U t ) -Π (q-1) f (U t ) ≤ Le 2 t + |e t ||f (U t ) -Π (q-1) f (U t )| ≤ (L + 1/2)e 2 t + 1 2 |f (U t ) -Π (q-1) f (U t )| 2 ,
where L is the Lipschitz constant on the second variable of g. By the Grönwall inequality, since e(0) = 0,

e 2 t ≤ t 0 |f (U s ) -Π (q-1) f (U s )| 2 e (2L+1)(t-s) ds
which yields an L 2 estimate of the error

∥u -U ∥ 2 L 2 ([0,T ],R n ) ≤ T 0 t 0 |f (U s ) -Π (q) f (U s )| 2 e (2L+1)(t-s) dsdt ≤ e (2 L +1)T 2L + 1 ∥f -Π (q-1) f ∥ 2 L 2 ([0,T ],R n ) (5.4)
In the following, we denote Π ζ a projector delivering an accuracy η ≥ 0 for the involved dynamics

∥f -Π (q-1) f ∥ L 2 ([0,T ],R n ) ≤ η.
From inequality (5.4), such a projector gives an accuracy in the solution U which is bounded by

∥u -U ∥ L 2 ([0,T ],R n ) ≤ e (2 L +1)T /2 (2L + 1) 1/2 η
The routines:

In the following, the routine

u η = solve_fwd(η; f, x)
yields an evolution u η = {u η t : t ∈ [0, T ]} that approximates the exact solution u of the forward dynamics (5.1) with a projector Π η . The routine

p η = solve_bwd(η; f, x)
works similarly for the backward dynamics ṗt = f (p t ) with p T = x.

The learning phase of the algorithm

The A-MSA algorithm consists in performing steps (4.3) and (4.4) of the original MSA but with a numerical time-integrators which give inexact trajectories. We carry the discussion assuming that we use the routines solve_fwd and solve_bwd introduced above. At every iteration k, both routines involves projectors Π k that deliver an accuracy ε k which is yet to be determined. The accuracy will be tightened as k increases in a way that still guarantees convergence to the exact, continuous solution (3.4) of the sampled PMP problem. We deduce its value later on from the convergence analysis of section 6.

Starting from an initial guess of the optimal control θ 0 , for each k ≥ 0 we first solve at accuracy ε k the forward dynamics

u θ k ,i,ε k = solve_fwd ε k ; (t, u) → f (u, θ k t ), x i , ∀i = 1, . . . , N.
Satisfying this accuracy may require to adapt the time discretization mesh and/or the numerical scheme as discussed in the previous section. An optimal control framework for adaptive neural ODEs

We then use the final state u θ k ,i,ε k T to compute an approximation of the backward dynamics at the same accuracy ε k for all i ∈ {1, . . . , N }

p θ k ,i,ε k = solve_bwd ε k ; (t, p) → -∇ u H(t, u θ k ,i,ε k t , p, θ k t ), -∇ u Φ(u θ k ,i,ε k T , y i ) .
This step may, again, require adaptivity to satisfy the target tolerance. To update the control, instead of finding the exact maximum like in (4.7), it is in fact sufficient to find controls that are away by a factor 0 < γ k ≤ 1 of the maximum and for which γ k → 1 as k → ∞. In other words, it suffices to find for all t ∈ [0, T ], a control

θ k+1 t ∈ Θ such that for all t ∈ [0, T ] 1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) ≥ γ k max θ∈Θ 1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ).
(5.5) Note that when γ k = 1, we fall back to the case where we look for the exact maximum. Algorithm 1 summarizes the whole procedure by describing a routine A-MSA[S, τ ] which computes a control dynamics θ that solves the extended PMP problem until the difference between two successive iterates is smaller than τ > 0. The algorithm also yields a numerical scheme associated to the parameter ε which is the accuracy of the forward and backward solvers used in the last iteration. These elements are then used for prediction as explained in the next section.

The practical implementation of the algorithm requires to specify the tolerances ε k and γ k . They can be theoretically derived from the convergence analysis of section 6. However, since the bounds of the analysis may not be sharp and that it involves quantities which are difficult to estimate in practice, devote section 8 to give indications on how to implement A-MSA in practice.

Predictions

The routine A-MSA[S, τ ] gives a discrete control θ and an accuracy ε. We also have a certain numerical scheme for the forward propagator, which can be interpreted as the final network architecture. As a result, for a given input data x, we can predict the output by first computing

û = solve_fwd ε; (t, u) → f (u, θ), x ,
and then taking

ŷ := g(û T )
as the approximation of the true y. In other words, our neural network is the mapping

N N : R n → R k Algorithm 1 Learning Algorithm: A-MSA[S N , τ ] → [ θ, ε] 1: Require: S N = {x i , y i } N i=1
, target tolerance τ 2: Set of internal parameters: θ 0 , ρ, maximum number of iterations k max . 3: k ← 0 4: J, J old ← J(θ 0 ) 5: ∆ ← τ + 1

u θ k ,i,ε k = solve_fwd ε k ; (t, u) → f (u, θ k t ), x i 11: p θ k ,i,ε k = solve_bwd ε k ; (t, p) → -∇xH(t, u θ k ,i,ε k t , p, θ k t ), -∇uΦ(u θ k ,i,ε k T , y i ) 12:
end for 13:

for t ∈ [0, T ] do ▷ In parallel 14:
Find θ k+1 t s.t.

1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) ≥ γ k max θ∈Θ 1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ) 15:
end for 16:

J old ← J; J ← J(θ k+1 ); ∆ ← J old -J 17: end while 18:

19: θ ← θ k+1 ; ε ← ε k 20: Output: [ θ, ε] x → N N (x) := g(û T ) = g solve_fwd ε; (t, u) → f (u, θ), x (t = T )
We summarize the prediction pipeline in Algorithm 2.

Algorithm 2 Prediction algorithm: N N (x) → ŷ 1: Input: -Observation x ∈ R n , -Neural network architecture [ θ, ε] (from A-MSA[S N , τ ]). 2: û ← solve_fwd ε; (t, u) → f (u, θ), x 3: Output: g(û T )
6 A priori convergence of A-MSA

In this section, we prove that, for suitably chosen tolerances (ε k ) k , the sequence of controls (θ k ) k≥0 computed with A-MSA is a minimizing sequence for problem (3.4), that is,

J S N (θ k ) k→∞ ----→ J * S N .
Our analysis is built upon the one presented in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF] for the exact E-MSA to which it is necessary to add perturbative arguments due to the inexact propagations.

We work with the same continuity assumptions as in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF] for the ideal case of exact propagations: (A1) Φ is twice continuously differentiable. Φ and ∇Φ satisfy the following Lipschitz condition: there exists K > 0 such that

|Φ(u) -Φ(u ′ )| + ∥∇Φ(u) -∇Φ(u ′ )∥ ≤ K∥u -u ′ ∥, ∀(u, u ′ ) ∈ R n .
(A2) For all t ∈ [0, T ] and θ ∈ Θ, u → f (u, θ) is twice continuously differentiable and satisfying the following Lipschitz condition: there exists

K ′ > 0 such that ∥f (u, θ)-f (u ′ , θ)∥+∥∇ u f (u, θ)-∇ u f (u ′ , θ)∥ 2 ≤ K∥u-u ′ ∥, ∀(u, u ′ ) ∈ R n ,
where ∥ • ∥ 2 denotes the induced 2-norm.

Intermediate problem:

We first study the following intermediate problem. Consider a fixed inputoutput pair (x, y). For a given control θ ∈ L ∞ ([0, T ], Θ), let u θ,x and p θ,x be the exact solutions of the forward and backward propagations using x and y in the initial and final conditions,

uθ,x t = f (u θ,x t , θ t ), u θ,x 0 = x, ṗθ,x t = -∇ u H(t, u θ,x t , p θ,x t , θ t ), p θ,x T = -∇ u Φ(u θ,x T , y),
Let u θ,x,ζ and p θ,x,ζ be the outputs of the same propagation but with accuracy ζ > 0,

u θ,x,ζ = solve_fwd (ζ; (t, u) → f (u, θ), x) p θ,x,ζ = solve_bwd ζ; (t, p) → -∇ u H(t, u θ,x,ζ t , p, θ t ), -∇ u Φ(u θ,x,ζ T , y) . (6.1)
Similarly, u φ,x,η and p φ,x,η denote the η-accurate solutions but for another control φ ∈ L ∞ ([0, T ], Θ).

The loss function associated to u θ,x,ζ is denoted

Loss(x, y, θ, ζ) := Φ(u θ,x,ζ T , y) + T 0 R(θ t )dt,
and similarly for u φ,x,η . In addition, we define the quantity

∆H x,ζ φ,θ (t) := H(t, u θ,ζ t , p θ,ζ t , φ t ) -H(t, u θ,ζ t , p θ,ζ t , θ t ), (6.2)
which will play an essential role in what follows. Lemma 6.1 is an important building-block in the proof of convergence of A-MSA. To not interrupt the flow of reading, we have deferred its proof to Appendix A. Lemma 6.1. Let θ, φ ∈ L ∞ ([0, T ], Θ) be two controls and let (x, y) ∈ R n ×R k be an input-output pair. Let (u θ,x,ζ , p θ,x,ζ ) and (u φ,x,η , p φ,x,η ) be solutions of the forward and backward associated problems with accuracy ζ and η like in (6.1). Then, there exist a constant C > 0 independent of θ, φ but dependent on T , ∥Π ζ ∥ and ∥Π η ∥ such that, if η ≤ 1,

Loss(x, y, φ, η) -Loss(x, y, θ, ζ) (6.3) ≤ - T 0 ∆H x,ζ φ,θ (t)dt + C (η + ζ) 2 + ∥∇ w (∆H x,ζ φ,θ )∥ 2 L 2 ([0,T ])

Convergence of A-MSA:

We can use Lemma 6.1 to prove our main convergence result which we give in the following theorem. It involves the quantity

λ 2 k := 1 N N i=1 ∥∇ w (∆H i,ε k θ k+1 ,θ k )∥ 2 L 2 ([0,T ]) = 1 N N i=1 T 0 ∥f (u θ k ,i,ε k t , θ k+1 t ) -f (u θ k ,i,ε k t , θ k t )∥ 2 dt + T 0 ∥∇ u H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) -∇ u H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k t )∥ 2 dt
Theorem 6.2. Let δ > 0 be a fixed constant and suppose that we run the A-MSA algorithm with the following accuracies:

• (ε k ) is a positive and decreasing sequence such that

0 ≤ ε k ≤ min ε k-1 , δλ 2 k 4C + ρN -1/2 λ k , ∀k ≥ 1.
• The maximization step (5.5) is performed with the parameter

γ k ≥ N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )/N δλ 2 k + N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )/N , ∀k ∈ N.
Then, if ρ > 2(C +2δ+ε 0 ), where C is the constant of Lemma 6.1, the sequence of controls (θ k ) of A-MSA is a minimizing sequence for problem (3.4), that is,

J S N (θ k ) → k→∞ J * S N .
Furthermore, λ 2 k → 0 as k → ∞ and we have

J S N (θ k+1 ) -J S N (θ k ) ≤ -κλ 2 k , κ := |C + 2δ - ρ 2 | > 0, ∀k ∈ N.
Proof We start by applying inequality (6.3) from Lemma 6.1 to the controls φ = θ k+1 and θ = θ k of the A-MSA algorithm, with accuracies η = ε k+1 and ζ = ε k and for the samples (x, y) = (x i , y i ). We then take the average over i and obtain

J S N (θ k+1 ) -J S N (θ k ) = 1 N N i=1 Loss(x i , y i , θ k+1 , ε k+1 ) -Loss(x i , y i , θ k , ε k ) ≤ - 1 N N i=1 T 0 ∆H i,ε k θ k+1 ,θ k (t)dt + C (ε k + ε k+1 ) 2 + λ 2 k (6.4)
We first do the proof when γ k = 1, ∀k ≥ 1. In this case, from the maximization step (5.5) in A-MSA,

1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t , uθ k ,i,ε k t , ṗθ k ,i,ε k t ) ≥ 1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k t , uθ k ,i,ε k t , ṗθ k ,i,ε k t ).
Recalling definition (4.6) for H, we reassemble the terms in the above inequality to derive

- 1 N N i=1 T 0 ∆H i,ε k θ k+1 ,θ k (t)dt ≤ - ρ 2N N i=1 T 0 ∥ uθ k ,i,ε k t -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 + ∥ ṗθ k ,i,ε k t -∇uH(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )| 2 + ρ 2N N i=1 T 0 ∥ uθ k ,i,ε k t -f (u θ k ,i,ε k t , θ k t )∥ 2 + ∥ ṗθ k ,i,ε k t -∇uH(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k t )∥ 2 dt (6.5)
Since by construction

uθ k ,i,ε k t = Πε k f (u θ k ,i,ε k t , θ k t ) and ∥f (u θ k ,i,ε k t , θ k t ) -Πε k f (u θ k ,i,ε k t , θ k t )∥ L 2 ([0,T ],R n ) ≤ ε k , we have by Cauchy-Schwartz inequality T 0 ∥ uθ k ,i,ε k t -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 dt = T 0 ∥f (u θ k ,i,ε k t , θ k t ) -Πf (u θ k ,i,ε k t , θ k t )∥ 2 + ∥f (u θ k ,i,ε k t , θ k t ) -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 -2 f (u θ k ,i,ε k t , θ k t ) -Πf (u θ k ,i,ε k t , θ k t ), f (u θ k ,i,ε k t , θ k t ) -f (u θ k ,i,ε k t , θ k+1 t ) dt ≥ T 0 ∥f (u θ k ,i,ε k t , θ k t ) -Πf (u θ k ,i,ε k t , θ k t )∥ 2 + ∥f (u θ k ,i,ε k t , θ k t ) -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 dt -2ε k T 0 ∥f (u θ k ,i,ε k t , θ k t ) -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 dt 1/2
Since a similar bound holds for

T 0 ∥ ṗθ k ,i,ε k t -∇uH(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )∥ 2 dt, we derive that - 1 N N i=1 T 0 ∆H i,ε k θ k+1 ,θ k (t)dt ≤ - ρ 2 λ 2 k + ρN -1/2 ε k λ k (6.6)
Using (6.6), we can further bound (6.4) to obtain

J S N (θ k+1 ) -J S N (θ k ) ≤ - ρ 2 λ 2 k + ρN -1/2 ε k λ k + C (ε k + ε k+1 ) 2 + λ 2 k ≤ C - ρ 2 λ 2 k + ε k (4C + ρN -1/2 λ k ) (6.7)
where we have used that (ε k ) is a decreasing sequence. Since, for a fixed δ > 0, we have chosen

ε k such that ε k ≤ δλ 2 k /(4C + ρN -1/2 λ k ), we infer that J S N (θ k+1 ) -J S N (θ k ) ≤ C + δ - ρ 2 λ 2 k . As a result, fixing ρ > 2(C + δ), J S N (θ k+1 ) -J S N (θ k ) ≤ -κλ 2 k , κ := C + δ - ρ 2 > 0.
Moreover, we can rearrange and sum over k = 0 to K the above expression to get

K k=0 λ 2 k ≤ κ-1 J S N (θ 0 ) -J S N (θ K+1 ) ≤ κ-1 J S N (θ 0 ) - inf θ∈L ∞ ([0,T ],Θ) J S N (θ) .
Therefore K k=0 λ 2 k < +∞ which implies that λ 2 k → k→∞ 0. Therefore, the A-MSA converges to a solution of the extended PMP.

Let us consider now the general case 0 < γ k ≤ 1. From the maximization step (5.5),

1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) ≥ γ k max θ∈Θ 1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ), ∀t ∈ [0, T ].
As a result, bound (6.5) receives an additional term and becomes,

- 1 N N i=1 T 0 ∆H i,ε k θ k+1 ,θ k (t)dt ≤ - ρ 2N N i=1 T 0 ∥ uθ k ,i,ε k t -f (u θ k ,i,ε k t , θ k+1 t )∥ 2 + ∥ ṗθ k ,i,ε k t -∇uH(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )∥ 2 dt + ρ 2N N i=1 T 0 ∥ uθ k ,i,ε k t -f (u θ k ,i,ε k t , θ k t )∥ 2 +∥ ṗθ k ,i,ε k t -∇uH(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k t )∥ 2 dt + γ -1 k -1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) .
It follows that bound (6.6) receives the same additional term

- 1 N N i=1 T 0 ∆H i,ε k θ k+1 ,θ k (t)dt ≤ - ρ 2 λ 2 k + ρN -1/2 ε k λ k + γ -1 k -1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) .
Therefore (6.7) becomes

J S N (θ k+1 ) -J S N (θ k ) ≤ C - ρ 2 λ 2 k +ε k (4C +ρN -1/2 λ k )+ γ -1 k -1 N N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t ) .
Choosing ε k like before and setting

γ k ≥ N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )/N δλ 2 k + N i=1 H(t, u θ k ,i,ε k t , p θ k ,i,ε k t , θ k+1 t )/N yields J S N (θ k+1 ) -J S N (θ k ) ≤ C + 2δ - ρ 2 λ 2 k .
We conclude along the same lines as before by fixing now ρ > 2(C +2δ), and inferring that

J S N (θ k+1 ) -J S N (θ k ) ≤ -κλ 2 k , κ := C + 2δ - ρ 2 > 0.
□ Before going to the next section, a few remarks on Theorem 6.2 are in order:

• Note that since the sequence of controls (θ k ) built by A-MSA is a minimizing sequence, it need not convergence in L ∞ ([0, T ], Θ) norm, nor does it necessarily contain a convergent subsequence. To guarantee that the minimizing sequence converges in norm, it is necessary to additional convexity and/or compactness properties. • The convergence result is a priori in the sense that it cannot be directly implemented for the following reasons:

-The bounds for (ε k ) and (γ k ) involve a constant C which is difficult to estimate in practice.

-The bounds may be suboptimal due to the construction of the proof.

-In practice, it is difficult to guarantee that the maximization step is performed within a certain fraction γ k of the actual maximum. • The main interest of the result lies in the fact that it reveals the necessity to tighten the accuracy in the time integration schemes and the maximization algorithms across the iterations in order to approach the continuous sampled optimal control problem (3.4).

Generalization error and convergence towards the fully continuous problem

The goal of this section is to connect the controls (θ k ) of the A-MSA Algorithm 1 with the PMP solutions of the fully continuous problem (3.2). For this, we rely on recent results from [START_REF] Han | A mean-field optimal control formulation of deep learning[END_REF] connecting the sampled PMP with the fully continuous PMP which we next briefly recall. We assume in the following that θ * and ϑ * S N are solutions of (3.2) and (3.4) such that the Hamiltonian step attains a maximum in the interior of Θ. Consequently, the continuous PMP solution θ * satisfies

F (θ * ) t := E µ ∇ θ H(t, u θ * t , p θ * t , θ * t ) = 0
for a.e. t ∈ [0, T ] and where F :

L ∞ ([0, T ], Θ) → L ∞ ([0, T ], R m ).
Similarly, an interior solution ϑ * S N of the sampled PMP is a random variable which satisfies

F S N (ϑ * S N ) t := 1 N N i=1 ∇ θ H(t, u ϑ * S N ,i t , p ϑ * S N ,i t , ϑ N,t ) = 0.
Under these assumptions, the theorem below, proven in [START_REF] Han | A mean-field optimal control formulation of deep learning[END_REF], describes the convergence of an interior solution ϑ * S N of the first order condition of the sampled PMP to an interior solution θ * of the continuous PMP. The additional local strong concavity assumption on the Hessian allows to guarantee that ϑ * S N is a global/local maximum of the sampled PMP. The result also guarantees convergence of loss function values.

Definition 7.1. For ρ > 0 and x ∈ L ∞ ([0, T ], Θ), define S ρ (x) := {y ∈ L ∞ ([0, T ], Θ) : ∥x -y∥ ≤ ρ}.
The mapping F is said to be stable on S ρ (x) if there exists a constant K ρ > 0 such that for all y, z ∈ S ρ (x),

∥y -z∥ L ∞ ([0,T ],Θ) ≤ K ρ ∥F (x) -F (y)∥ L ∞ ([0,T ],R m ) .
Theorem 7.2 (Theorem 6, Corollary 1 and 2 from [START_REF] Han | A mean-field optimal control formulation of deep learning[END_REF]). Let F be a mapping which is stable on S ρ (θ * ) for some ρ > 0 and θ * a solution of F = 0. Then there exists positive constants s 0 , C, K 1 , K 2 and ρ 1 < ρ and a random variable

ϑ * S N ∈ S ρ1 (θ * ) such that P[∥θ * -ϑ * S N ∥ ≥ Cs] ≤ 4 exp - N s 2 K 1 + K 2 s
, s ∈ (0, s 0 ] and ϑ * S N → θ * as N → ∞ in probability. Moreover, there exists constants

K ′ 1 , K ′ 2 such that P[ J (θ * ) -J (ϑ * S N ) ≥ s] ≤ 4 exp - N s 2 K ′ 1 + K ′ 2 s
, s ∈ (0, s 0 ] (7.1)

In addition, if for all t ∈ [0, T ], the Hessian E µ ∇ 2 θ,θ H(u θ * t , p θ * t , θ * t ) + λ 0 I is negative definite, then ϑ * S N is also a strict local maximum of the sampled Hamiltonian ϑ → 1 N N i=1 H(u ϑ * S N ,i t , p ϑ * S N
,i t , ϑ). In particular, if the sampled Hamiltonian has a unique maximizer, then ϑ * S N is a solution of the sampled PMP (4.2) with the same high probability.

Corollary 7.3. If the sampled Hamiltonian has a unique maximizer ϑ * S N , then there exists constants K1 and K2 such that for any ε > 0, there exists k(ε) ∈ N such that for k ≥ k(ε), the control θ k at iteration k of A-MSA satisfies

P |J S N (θ k ) -J (θ * )| ≥ 3ε ≤ 4 exp - N ε 2 K1 + K2 ε Proof We use that P |J S N (θ k ) -J (θ * )| ≥ 3ε ≤ P |J S N (θ k ) -J S N (ϑ * S N )| ≥ ε + P |J S N (ϑ * S N ) -J (ϑ * S N )| ≥ ε + P |J (ϑ * S N ) -J (θ * )| ≥ ε ,
and then bound each term as explained next. Let ε > 0 be fixed. With high probability, ϑ * S N is a solution to the sampled PMP (3.4), that is

J S N (ϑ * S N ) = J * S N . By Theorem 6.2, there exists k(ε) ∈ N such that for k ≥ k(ε), the control θ k at iteration k of A-MSA satisfies |J S N (θ k ) -J S N (ϑ * S N )| ≤ ε,
which bounds the first term. For the second term, by the infinite dimensional Hoeffding's inequality (see [51,

Corol- lary 2]), |J S N (ϑ * S N ) -J (ϑ * S N )| ≤ ε with probability 1 -exp(-N ε 2 /( K′ 1 + K′ 2 ε))
. The third term is also bounded by ε with probability 4 exp -N ε 2 K ′ 1 +K ′ 2 ε thanks to (7.1). □

Guidelines for numerical implementation

Algorithm 1 should in principle be implemented with tolerances ε k and γ k given in Theorem 6.2 to ensure convergence. These tolerances are in practice hard to guarantee, and the convergence analysis that leads to these quantities may be suboptimal. A trade-off between the theory and numerical implementation consists simply in replacing the forward and backward propagations at prescribed accuracies by propagations where the time steps are progressively refined across the iterations. It is therefore necessary to prescribe a refinement strategy. For the Hamiltonian maximization step (line 13 of Algorithm 1), the best that one can do is to use an efficient optimizer of nonconvex problems.

In addition to these considerations, note that Algorithm 1 can be parallelized at several parts: for each sample i, the forward and backward propagations can be performed in parallel (see the for loop of line 8). Each propagation can additionally be parallelized by means of parallel in time algorithms as in [START_REF] Günther | Layerparallel training of deep residual neural networks[END_REF] (see [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF][START_REF] Maday | An adaptive parareal algorithm[END_REF] for some selected references on the topic). In addition, the Hamiltonian maximization step can also be parallelized (see for loop of line 12 of Algorithm 1). For our development, we have chosen to parallelize the Maximization step, since it is way more expensive in terms of computational time than the propagation step.

Numerical experiments

In this section, we present some numerical experiments, aiming primarily at illustrating the behavior of the A-MSA algorithm. We give special focus on studying the performance in the learning phase in terms of the value of the loss function. We also discuss computing times and examine expressivity in terms of generalization errors. Note that the final quality of approximation depends on the number of samples, on the network depth, and on the optimizers and their starting guesses. We proceed by increasing levels of difficulty in the numerical examples in order to disentangle the effect of each of the different aspects (note however that their effects are usually mixed in non-trivial applications). The implementation has been done with Python 3 and NumPy. Particular attention has been paid to vectorize most of the operations in order to delegate as much as possible the looping to internal, highly optimized C and Fortran functions. We have observed that the time to perform the ODE propagations is negligible with respect to the Hamiltonian maximization (line 12 of Algorithm 1). We have thus parallelized this step with MPI.

Approximation of the sine function

In this example, inspired from [31, Section 6], the task is to approximate the graph of the sine function,

y : X = [-π, π] → Y = [-1, 1],
x → y(x) = sin(x).

Note that here the dimensions n and k of the domain and label sets are n = k = 1. The continuous learning problem (see (3.2)) is to minimize the L 2 (X)

approximation error over the controls θ ∈ L ∞ ([0, T ], Θ), namely inf θ∈L ∞ ([0,T ],Θ) 1 2 X |y(x) -g(u θ,x T )| 2 dx.
For every point x ∈ X, u θ,x T is the final time state of the ODE

uθ,x t = tanh(A t • u θ,x t + b t ), ∀t ∈ (0, T ] u θ,x 0 = x(1, . . . , 1) T ∈ R d .
Here,

θ t = (A t , b t ) ∈ R d×d × R d ∼ R d 2 +d are the parameters to optimize. We constraint them to lie in [α min , α max ] = [-1, 1] so Θ = [-1, 1] d 2 +d
. The function g is defined as

g : R d → R, z = (z 1 , . . . , d ) → g(z) := 1 d d i=1 z i ,
and note that we have chosen µ = dx as the Lebesgue measure, f = tanh as the activation function, Φ(u θ,x T , y) = 1 2 |y(x) -g(u θ,x T )| 2 and we do not have any regularisation (R = 0).

In the sampled version which we consider in our experiments, we work with a set S N = {(x i , y i )} N i=1 of pairs generated with equidistant data points

x i = -π + 2π i -1 N -1 , y i = sin(x i ), i ∈ {1, . . . , N },
and perform empirical risk minimization taking the uniform distribution µ N = 1 N N i=1 δ (xi,yi) (see (3.4)). We solve this problem with the above described A-MSA algorithm where the final time is set to T = 5, the penalization parameter is set to ρ = 5, and the number of neurons per layer to d = 3. We use a basic explicit Euler scheme as a time integrator, which induces a ResNet architecture. We progressively refine the time discretization from L = 3 to L = 32 time steps (in other words, we increase the network depth from 3 to 32 layers) using three different strategies: (A1) Abrupt refinement: We refine from shallow (L = 3) to deep (L = 32) at iteration k = 250. (A2) Fast refinement: We add 10 layers every 50 iterations. (A3) Slow refinement: We add 10 layers every 100 iterations. Our goal is to compare the behavior of these refinement strategies with the non-adaptive training of a shallow and a deep network having L = 3 and L = 32 layers respectively.

At each iteration k ≥ 0 of A-MSA, the maximization step is performed with NumPy's L-BFGS-B optimizer, which is an iterative algorithm suited for nonconvex optimization problems. The optimizer is initialized by picking the best candidate among a list Θ 0 of vectors, including the last best control θ k best = arg min 0≤i<k J S N (θ i ) among all previous iterations before k and random perturbations of it, precisely

Θ 0 := 5 q=0 25 i=1 θ k best + 10 -2q rand(α min , α max ), 10 -2q rand(α min , α max )
where rand(α min , α max ) is a uniform real distribution function in the interval [α min , α max ]. As a result, the final output of the algorithm is random due to the choice of the initial guess. For this reason, we make r = 20 runs of A-MSA, each one consisting in k max = 800 iterations.

In Figure 1, we first fix N = 20 samples and we plot the convergence history of the train loss J (run=r) S N (θ k ) (see red colors). We also plot the generalization error via the test loss across the iterations k (see blue colors). The test loss is the quantity J S N (θ k ), defined in (3.4), where we use a test set S N different from the training set S N . S N is taken uniformly random at each iteration 1 ≤ k ≤ k max . Since the history depends on each run r, the figure shows some statistics related to the repetitions. The continuous red curve shows the average value J (av)

S N ,k := (1/r) r r=1 J (run=r) S N
(θ k ) over the runs and the red diffuse color shows its distribution around the average. The same presentation is given in blue for the generalization errors. We can first observe that the convergence history is, in average, relatively similar for all architectures. The average value of the loss function does not go below 10 -2 for the shallow network (figure 1a), and the value reached by the deep network is only slightly better (figure 1b): at the end of the iterations, the cost function is in average below 7.10 -2 . Similar values are reached by any of our three adaptive strategies (see figures 1d, 1e and 1f). As a consequence, if we take this average value as the performance indicator, we may conclude that the approximation power of the deep network is only marginally better than the shallow one. We may also conclude that the adaptive training strategy reduces to some extend the computational time for training (as we illustrate further on in section 9.4) but it does not bring extra approximation power. However, note that this indicator is not entirely appropriate because it does not correspond to any realized convergence history and, more importantly, because it does not inform about the best performance that we have at hand. It is crucial to remark that the variance around the average convergence value is significantly larger in the adaptive strategies compared to the non adaptive ones. This is particularly true for (A2) and (A3). Figure 1c shows that the minimal value min 1≤k≤kmax J (run=r) S N (θ k ) reached during the k max = 800 iterations of each run is, in average, one order of magnitude lower than the one of the non-adaptive coarse and the non-adaptive deep network. This result illustrates that the adaptive strategy significantly contributes to reach better minima, thus allowing to benefit better from the higher approximation power of the deeper neural network. (θ k ) for every r. We see that the best over all realizations yields a very satisfactory approximation. Statistics on the generalization errors of these best controls are given in Figure 2c We can next examine in Figure 3 the impact of the number N of data samples in the generalization errors of the best controls. As one can expect, the generalization errors decreases when the number N of samples increases. For N ≥ 20, the errors with the shallow network stagnate before reaching 10 -3 in average regardless of the number of samples. This illustrates the limitations in the approximation power of the shallow network in the current example. The non-adaptive deep network (L = 32) presents a behavior which is only better for a reduced number of samples N ≤ 8 but overall its performance is relatively similar to the coarse neural network. This comes from the difficulty in finding good quality optimizers in networks involving many coefficients. If we now adaptively train the deep network with strategy (A2), the generalization errors are in general lower than the ones reached by the non-adaptive strategies, especially when N ≥ 15. The result illustrates that the adaptive strategy gives in average a better performance in terms of approximation. This is obtained at a reduced computing time as we show further on.

Noisy step function

We next examine the robustness of the method against noise.We consider the graph of the step function, and we seek to approximate it when the given data are noisy. For each sample i, we get the pair (x i , y i ) with

y : X = [-1, 1] → Y = [-0.5, 0.5], x → y(x) = 0.5 if x ≤ 0 -0.5 if x > 0
y i = y(x i ) + ε i , i = 1, . . . , N,
and the ε i follow a uniform distribution in [-0.2, 0.2].

We consider the same setting as before (T = 5, ρ = 5, d = 3) and fix N = 800. Figure 4 shows the convergence history of the training with the shallow, deep and adaptively deep neural networks. We consider the same adaptive strategies as before. All approaches reach more or less the same minimal value in the training (see Figure 4c. However, we observe that (A2) and (A3) give better generalization errors than the non-adaptive deep architecture (see Figure 5c). Note however that this superiority is not very large and all methods produce excellent reconstruction results as Figure 5 illustrates. We think that this is due to the fact that in this case all methods find configurations that are very close the global optimum, which is here given by the approximation error of the mapping x → y(x).

Classification

As a last example, we consider a simple 2D classification problem where the function to approximate is y :

X = [-1, 1] 2 → Y = {0, 1} s.t. y(x 1 , x 2 ) = 1 if x 2 1 + x 2 2 ≤ (0.5) 2 , 0 otherwise.
The setting is the same as in the above examples but with slight changes. Each layer has now d = 2 × 3 = 6 neurons per layers, the constraint on the controls is set to [α min , α max ] = [-2, 2] and the output function g is defined as A-MSA yields the same quality of approximation as the non-adaptive deep neural network.

g : R d → R, z = (z 1 , . . . , z d ) → g(z) := H 1 d d i=1 z i

Computational times

We finish the section on numerical tests by examining the computational time in the learning approach. We present runtimes only for the noisy step example for the sake of brevity (similar results are obtained for the other tests). We consider three different criteria:

• Complexity index (see Figure 8a): we estimate the number of operations by measuring the runtime of a sequential run. • Parallel runtime index (see Figure 8b): we measure the runtime of running in parallel the Hamiltonian maximization step of the learning algorithm. We place ourselves in a scenario with no constraints in the computing ressources and allocate one processor per layer. Note that the number of layers is increased in the adaptive approach so we use more and more ressources as the algorithm makes refinements. • Energy consumption index (see Figure 8a): we estimate the total energy consumption of a parallel run as follows. For each iteration of the algorithm, we add the corresponding computing times of all the processors involved in that iteration. This gives an estimate in cpu.seconds of the energy consumption for each iteration. By summing over all iterations, we obtain the final estimate. This quantity is in fact equal to the cumulative sequential runtime which is given on Figure 8a. Our runtimes have been measured on a cluster and we have observed that its occupation greatly impacts on the runtime results. To illustrate this issue, we show runtimes for the training of the shallow neural network in an empty node and a busy node, i.e. a node where others jobs are running in parallel. We have run the rest of the examples on empty nodes.

If we first consider Figure 8a, we see that the adaptive training strategy performs better than the non-adaptive one regarding the sequential runtimes and the energy consumption. This is comes as no surprise since the adaptive strategy performs less operations at the beginning of the iterations since they involve less layers. When computations are run in parallel, Figure 8b) shows that the run times are all very similar for all strategies. This is due to the fact that the tasks are well-balanced between processors since each of them does computations for one layer. We also observe that runtimes are greatly affected by the occupation of the cluster.

Conclusion

We have proven that the A-MSA algorithm converges to an underlying continuous learning problem in the limit of the time discretization and of the the number of samples. The convergence analysis requires that the time propagations and Hamiltonian maximizations are performed at each step with increasingly tight accuracy tolerances. Since a sharp estimation of these tolerances is difficult to obtain, we have implemented A-MSA with certain refinement strategies in order to examine its potential. The numerical experiments reveal that the adaptive strategy helps to benefit in practice from the higher approximation properties of deep networks by mitigating overparametrization issues: our results show that adaptivity increases the chances to find better quality minimizers compared to the non-adaptive training of deep neural networks. In addition, it appears that the adaptive strategy is clearly more performant in terms of complexity and energy consumption compared to the non-adaptive training of a deep neural network. Some limitations and natural questions that arise for future works are the following:

• Our theoretical results involve numerous bounds and some of them might be suboptimal. Certain constants involved are also difficult to estimate in practice. As a result, the theoretical convergence result cannot be used to derive fully explicit adaptive strategies for implementation. However it gives hints on how to proceed in practice, and we propose practical guidelines to perform adaptivity. In future works, it would be desirable to find a sharper analysis with constants that are easier to compute in order to obtain a fully closed pipeline. This may however require very different arguments as our current ones. • As we will explain in the next section, the optimal control point of view only describes certain classes of neural networks since it introduces certain restrictions in the architecture. As a result, our adaptive setting cannot automate the training of all existing classes of neural networks. So a natural question is how to extend the current ideas to other training techniques that are not based on optimal control. • One important question is whether the flow map of a given dynamical system can represent the data at hand. The same issue also arises in the classical approach of deep learning in terms of the optimality of the choice of the network. From the optimal control perspective, this representability problem can be viewed as a controllability problem. It would be interesting to transfer controllability results to the current deep learning framework to address this question.

Proposition A.1. We have the following bounds for all t ∈ [0, T ],

∥p θ,ζ t ∥ ≤ K ′ := Ke ∥Π ζ ∥KT (A.1) ∥δu t ∥ ≤ K δu := √ T (η + ζ) + T 0 ∥f (s, u θ,ζ s , φ s ) -f (s, u θ,ζ s , θ s )∥ds e KT (A.2) ∥δp t ∥ ≤ K δp := √ T e 2KT (1 + K(1 + K ′ T ))(η + ζ) + e 2KT K(1 + K ′ T ) T 0 ∥f (s, u θ,ζ s , φ s ) -f (s, u θ,ζ s , θ s )∥ds + e KT T 0 ∥∇ u H(s, u θ,ζ s , p θ,ζ s , φ s ) -∇ u H(s, u θ,η s , p θ,η s , θ s )∥ds (A.3)
Proof of Proposition A.1 Inequality (A.1) follows from the fact that

ṗθ,ζ t = Π ζ -p θ,ζ t • ∇uf (u θ,ζ t , θ t )
thus taking the scalar product with p θ,ζ t and using the Cauchy-Schwarz inequality, we derive the evolution 1 2

d dt ∥p θ,ζ t ∥ 2 ≤ ∥Π ζ ∥∥∇uf (u θ,ζ t , θ t )∥ 2 ∥p θ,ζ t ∥ 2 , ∥p θ,ζ T ∥ 2 = ∥∇uΦ(u θ,ζ T , y)∥ 2

6 : 7 : 8 :ε k and γ k 9 :

 6789 while ∆ > τ or k ≤ k max do Update for i in {1, . . . , N } do ▷ In parallel 10:

Fig. 1 :

 1 Fig. 1: Sine function. Training loss J (run=r) S N (θ k ) (red) and test loss

Figure 2

 2 Figure 2 shows the reconstruction of the sine function with the controls performing best in each training run, namely, realizing min 1≤k≤kmax J (run=r) S N

Fig. 2 :

 2 Fig. 2: Sine function: Predictions of r = 20 runs and average prediction.

Fig. 3 :

 3 Fig. 3: Sine case -Generalization errors vs N for coarse, adaptive (A2) and thin networks

Fig. 4 :

 4 Fig. 4: Noisy step: Loss J (run=r) N (θ k ) for N = 800 and r = 20 runs.

Fig. 5 :

 5 Fig. 5: Noisy step: Predictions of some runs.

  where H is the usual Heaviside step function. The figure layout of the results is the same as in the other examples. Figures6 and 7show statistics on the loss function and examples of reconstruction. The blue-to-red scalar field in Fig5is the output of the trained network composed with a treshold filter F : x → H(x -0.5) from 1024 inputs points. The black crosses (×) and white dots (•) are the N = 800 training points used for each case, the first refers to points with z = 0, the second to z = 1. Essentially, we observe that

Fig. 6 :

 6 Fig. 6: Classification: Loss J (run=r) N (θ k ) for N = 800 and 20 runs.

Fig. 7 :

 7 Fig. 7: Classification: Predictions of r = 20 runs and average prediction. Data labels: × : 0 , • : 1 .

  (a) Complexity and energy consumption. (b) Parallel runtime.

Fig. 8 :

 8 Fig. 8: Computational times for the noisy step example.

1 .

 1 and the result follows by applying the Grönwall's inequality and the fact that ∥∇uΦ(u θ,ζ T , y)∥ ≤ K and ∥∇uf (u θ,ζ t , θ t )∥ 2 ≤ K by hypothesis (A1) and (A2). To derive inequality (A.2), we start fromδ ut = uφ,η t -uθ,ζ t = Πηf (u φ,η t , φ t ) -Π ζ f (u θ,ζ t , θ t ).Integrating between 0 and t ≤ T ,δu t = t 0 Πηf (s, u φ,η s , φs) -Π ζ f (s, u θ,ζ s , θs)dsTaking norms,∥δu t ∥ ≤ t 0 ∥f (s, u φ,η s , φs) -Πηf (s, u φ,η s , φs)∥ds + t 0 ∥f (s, u θ,ζ s , θs) -Π ζ f (s, u θ,ζ s , θs)∥ds + t 0 ∥f (s, u φ,η s , φs) -f (s, u θ,ζ s , θs)∥ds ≤ √ T (η + ζ) + T 0 ∥f (s, u φ,η s , φs) -f (s, u θ,ζ s , φs)∥ds + T 0 ∥f (s, u θ,ζ s , φs) -f (s, u θ,ζ s , θs)∥ds u θ,ζ s , φs) -f (s, u θ,ζ s , θs)∥dsThus by Grönwall's inequality,∥δu t ∥ ≤ K δu := √ T (η + ζ) + T 0 ∥f (s, u θ,ζ s , φs) -f (s, u θ,ζ s , θs)∥ds e KT Proceeding in a similar manner with δ ṗt , by integrating between T and t, we derive t T δ ṗsds = δp t -δp T = t T Πη∇uH(s, u φ,η s , p φ,η s , φs) -Π ζ ∇uH(s, u θ,ζ s , p θ,ζ s , θs)ds and similarly as before, we infer that ∥δp t ∥ ≤ √ T (η + ζ) + K∥δu T ∥ + KK ′ , u θ,ζ s , p θ,ζ s , φs) -∇uH(s, u θ,ζ s , p θ,ζ s , θs)∥ds Using bounds (A.2) and (A.3), and applying the Gronwall inequality, we derive the announced estimate∥δp t ∥ ≤ K δp with K δp = e KT √ T (η + ζ) + KK δu (1 + K ′ T ) + T 0 ∥∇uH(s, u θ,ζ s , p θ,ζ s , φs) -∇uH(s, u θ,ζ s , p θ,ζ s , θs)∥ds = √ T e 2KT (1 + K(1 + K ′ T ))(η + ζ) + e 2KT K(1 + K ′ T ) T 0 ∥f (s, u θ,ζ s , φs) -f (s, u θ,ζ s , θs)∥ds + e KT T 0 ∥∇uH(s, u θ,ζ s , p θ,ζ s , φs) -∇uH(s, u θ,ζ s , p θ,ζ s , θs)∥ds □ We can now prove Lemma 6.Proof of Lemma 6.1 By definition (4.1) of the Hamiltonian, for any θ ∈ L ∞ ([0, T ], Θ), I(u θ , p θ , θ) := T 0 p θ t • f (u θ t , θ t ) -H(t, u θ t , p θ t , θ t ) -L(θ t ) dt = 0.Thus 0 = I(u φ,η , p φ,η , φ) -I(u θ,ζ , p θ,ζ , θ)

  φ,η t , φ t ) -p θ,ζ (u φ,η t , φ t ) -Πηf (u φ,η t , φ t ) and similarly for e f,θ,ζ By integration by parts, and the fact that ṗθ,ζ t = -Π ζ ∇uH(t, u θ,ζ t , p θ,ζ t , θ t ), we have I 1,1,1 = p θ,ζ Π ζ ∇wH(t, w θ,ζ t , θ t ) • δw t dt,where we have used the shorthand notation w = (u, p) in the last line. This notation will also be used in what follows.We next address the intergral I 1,1,2 . By integration by parts,Πη∇wH(t, w φ,η t , φ t ) -Π ζ ∇wH(t,w θ,ζ t , θ t ) • δw t dt Adding and subtracting ∇wH(t, w φ,η t , φ t ) and ∇wH(t, w θ,ζ t , θ t ) to the last formula, By applying Taylor's theorem around w = w θ,ζ = (u θ,ζ t , p θ,ζ t ) to the function w → ∇wH(t, w, φ t ), there exists r 1 (t) ∈ [0, 1] such that ∇wH(t, w φ,η t , φ t ) = ∇wH(t, w θ,ζ t , φ t ) + δw t • ∇ 2 w H(t, w θ,ζ t + r 1 (t)δw t , φ t ) (t)δw t , φ t ) • δw t dtSince δu 0 = 0, the boundary terms from I 1,1,1 and I 1,1,2 are there existsr 2 ∈ [0, 1] such that ∇Φ(u φ,η T ) = ∇Φ(u θ,ζ T ) + δu t • ∇ 2 Φ(u θ,ζ T + r 2 δu T ). δu t • ∇ 2 Φ(u θ,ζ T + r 3 δu t )• δu T , we finally get the expression for the the boundary terms from I 1,1,1 andI 1,1,2 • ∇ 2 Φ(u θ,ζ T + r 2 δu T ) -∇ 2 Φ(u θ,ζ T + r 3 δu η T ) • δu TWe now turn to I 2 . By Taylor's theorem, there exists r 4 (t) ∈ [0, 1] such that

	Denoting
		e f,φ,η
	t	, we have
	t • Π ζ f (u θ,ζ t , θ t )dt • f (u θ,ζ t , θ t )dt + T 0 p φ,η t Πη∇wH(t, w φ,η t , φ t ) -p θ,ζ • Πηf (u φ,η t :=I 1,1 δp t • δu t p φ,η t I 1,1,2 = I 1 = T 0 Thus 1 2 T 0 + 1 2 T t , φ t ) -∇wH(t, w φ,η • e f,φ,η t t , φ t ) • δw t dt -p θ,ζ t • e f,θ,ζ t 0 :=I 1,2 In addition, since uφ,η t = Πηf (u φ,η + 1 2 T ∇wH(t, w θ,ζ t , θ t ) -Π ζ ∇wH(t, w θ,ζ t , θ t ) • δw t dt 0 t , φ t ) and similarly for θ, it follows that I 1,1 = T 0 p φ,η t • uφ,η t -p θ,ζ t • uθ,ζ t dt = T 0 p θ,ζ t • δ ut + δp t • uθ,ζ t dt :=I 1,1,1 + T 0 + 1 2 T ∇wH(t, w θ,ζ t , φ t ) -∇wH(t, w θ,ζ t , θ t ) • δw t dt 0 δp t • δ ut dt :=I 1,1,2 where + 1 2 T 0 δw t • ∇ 2 w H(t, w θ,ζ t + r 1 p θ,ζ t + 1 2 δp t • δu t T 0 = p θ,ζ T + 1 2 δp T • δu T = -∇Φ(u θ,ζ T ) • δu T -1 2 ∇Φ(u φ,η Therefore δu = p θ,ζ t • δu t T 0 + T 0 I 1,1,2 = 1 2 T 0 δp t • δ ut dt + 1 2 T 0 δp t • δ ut dt = 1 2 δp t • δu t T 0 -1 2 T 0 δ ṗt • δu t dt + 1 T p θ,ζ t + 1 2 δp t • δu t T 0 = -∇Φ(u θ,ζ T ) • δu T -1 2 δu p θ,ζ t + 1 2 δp t • δu t T 0 = Φ(u θ,ζ T ) -Φ(u φ,η T )-δp t • δ ut dt 2 0 = 1 2 δp t • δu t T 0 + 1 2 T 0 1 2 δu T I 2 = T H(t, w φ,η t , φ t ) -H(t, w θ,ζ t , θ t )dt 0 :=I 1 -T 0 H(t, u φ,η t , p φ,η t , φ t ) -H(t, u θ,ζ t , p θ,ζ t , θ t )dt :=I 2 I 1,1,2 = 1 2 δp t • δu t T 0 + 1 2 T Πη∇wH(t, w φ,η t , φ t ) -∇wH(t, w φ,η t , φ t ) • δw t dt = T 0 ∆H x,ζ φ,θ (t)dt + T ∇wH(t, w θ,ζ t , φ t ) • δw t dt 0 0 + 1 2 1 + 2	dt
		:=I 3

-T 0 R(φ t ) -R(θ t )dt t := f t := u φ,η t -u θ,ζ

t , and δp t := p φ,η

t -p θ,ζ t . t • δu t T 0 + T 0 δp t • uθ,ζ t -ṗθ,ζ t • δu t dt = p θ,ζ t • δu t T 0 + T 0 δu t • Π ζ ∇uH(t, u θ,ζ t , p θ,ζ t , θ t ) + δp t • Π ζ ∇pH(t, u θ,ζ t , p θ,ζ t , θ t )dt T 0 ∇wH(t, w θ,ζ t , θ t ) -Π ζ ∇wH(t, w θ,ζ t , θ t ) • δw t dt + 1 2 T 0 ∇wH(t, w φ,η t , φ t ) -∇wH(t, w θ,ζ t , θ t ) • δw t dt T ) -∇Φ(u θ,ζ T ) • δu T By applying Taylor's theorem around z = u θ,ζ t to the function z → ∇Φ(z), T • ∇ 2 Φ(u θ,ζ T + r 2 δu T ) • δu T

Since, again by Taylor's theorem, there exists r 3 ∈ [0, 1] such that

Φ(p φ,η T ) = Φ(p θ,ζ T ) + δu t • ∇Φ(u θ,ζ T ) + 1 2 T 0 δw t • ∇

2 w H(t, w θ,ζ t + r 4 (t)δw t , φ t ) • δw t dt, where we remind that, as defined in equation (6.2), ∆H x,ζ φ,θ (t) := H(t, w θ,ζ t , φ t ) -H(t, w θ,ζ t , θ t ). Injecting the above expressions for the terms of I 1 and I 2 into equation

In fact, ResNets were originally defined for h = 1. The "augmented" definition with the paramater h was introduced in[START_REF] Haber | Stable architectures for deep neural networks[END_REF].

Acknowledgments and Disclosure of Funding

This research has been funded by the Emergences Project of the Paris city council called "Models and Measures".

Declarations

Conflict of interest The authors declare no competing interests

A Proof of Lemma 6.1

The proof of Lemma 6.1 makes use of certain bounds which we gather in the following Proposition A.1. They are derived by an immediate extension of the proofs given in [START_REF] Li | Maximum principle based algorithms for deep learning[END_REF] so we just sketch the proof for self-completness of the current work. An optimal control framework for adaptive neural ODEs (A.4), passing the terms corresponding to the loss function to the left hand side, and adding and subtracting

We next derive a bound for the difference Loss(x, y, φ) -Loss(x, y, θ). We proceed to bound all the terms by order of appearance in the above formula. We will sometimes use that by (A.2), (A.3) and Jensen's inequality,

which yields

Term (T1)+(T5):

By the Cauchy-Schwartz inequality and (A.6), we can bound term (T1) by

Term (T2):

Using (A.1) and the fact that ∥e

Terms (T3) and (T8):

By the Cauchy-Schwarz inequality, we can bound term (T3)

where we have used that the propagator Πη guarantees

Term (T4): Term (T6), (T9) and (T10):

By assumptions (A1), (A2), all second derivative terms are bounded element-wise by some constant K. Hence, we have ∥δw t • A • δw t ∥ ≲ ∥δw t ∥ 2 for A being a second derivative matrix.

Summary:

By using the above inequalities, in formula (Diff-Loss), we deduce that there exists a constant C > 0 depending on T , ∥Π ζ ∥ and ∥Πη∥ such that, if η ≤ 1, Loss(x, y, φ) -Loss(x, y, θ) ≤ -

Note that the term (η + ζ) 2 accounts for the fact that we are not solving exactly the dynamics. □