Camille Castera
email: camille.castera@protonmail.com

Jérôme Bolte

Cédric Févotte

Edouard Pauwels

An Inertial Newton Algorithm for Deep Learning

Keywords: deep learning, non-convex optimization, second-order methods, dynamical systems, stochastic optimization

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Can we devise a learning algorithm for general deep neural networks (DNNs) featuring inertia and Newtonian directional intelligence only by means of backpropagation? In an optimization jargon: can we use second-order ideas in time and space for non-smooth non-convex optimization by only using a subgradient oracle? Before providing some answers to this question, let us have a glimpse at some fundamental optimization algorithms for training DNNs.

The backpropagation algorithm is, to this day, the fundamental block to compute gradients in deep learning (DL). It is used in most instances of the Stochastic Gradient Descent (SGD) algorithm [START_REF] Robbins | A stochastic approximation method[END_REF]. The latter is powerful, flexible, capable of handling large-size problems, noise, and further comes with theoretical guarantees of many kinds. We refer to [START_REF] Bottou | The tradeoffs of large scale learning[END_REF], [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF] in a convex machine learning context and [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF] for a recent account highlighting the importance of DL applications and their challenges. In the non-convex setting, recent works of [START_REF] Adil | Opérateurs monotones aléatoires et application à l'optimisation stochastique[END_REF], [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF] follow the Ordinary Differential Equations (ODE) approach introduced in [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF], and further developed in [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF], [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF], [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF], [START_REF] Borkar | Stochastic approximation: A dynamical systems viewpoint[END_REF]. Two research directions have been explored in order to improve SGD's training efficiency:

-using local geometry of empirical loss functions to obtain steeper descent directions, -using past steps history to design larger step-sizes in the present.

The first approach is akin to quasi-Newton methods while the second revolves around Polyak's inertial method [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF]. The latter is inspired by the following appealing mechanical thought-experiment. Consider a heavy ball evolving on the graph of the loss function (the loss function's landscape), subject to gravity and stabilized by some friction effects. Friction generates energy dissipation, so that the particle will eventually reach a steady state which one hopes to be a local minimum. These two approaches are already present in the DL literature: among the most popular algorithms for training DNNs, ADAGRAD [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF] features local geometrical aspects while ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] combines inertial ideas with stepsizes similar to the ones of ADAGRAD. Stochastic Newton and quasi-Newton algorithms have been considered by [START_REF] Martens | Deep learning via Hessian-free optimization[END_REF], [START_REF] Byrd | On the use of stochastic Hessian information in optimization methods for machine learning[END_REF][START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF] and recently reported performing efficiently on several problems [START_REF] Berahas | An investigation of Newton-sketch and subsampled Newton methods[END_REF][START_REF] Xu | Second-order optimization for non-convex machine learning: an empirical study[END_REF]. The work of [START_REF] Wilson | The marginal value of adaptive gradient methods in machine learning[END_REF] demonstrates that carefully tuned SGD and heavy-ball algorithms are competitive with concurrent methods.

However, deviating from the simplicity of SGD also comes with major challenges because of the high dimensionality of DL problems and the severe absence of regularity in DL (differential regularity is generally absent, but even weaker regularity such as semi-convexity or Clarke regularity are not always available). All sorts of practical and theoretical hardships are met: computing and even defining the Hessian is delicate, inverting them is unthinkable to this day, first and second-order Taylor approximations are unavailable due to non-smoothness, and finally one has to deal with "shocks" which are inherent to inertial approaches in a non-smooth context ("corners" and "walls" in the landscape of the loss function generate velocity discontinuity). This makes in particular the study of the popular algorithms ADAGRAD and ADAM in full generality quite difficult, with recent progresses reported in [START_REF] Barakat | Convergence and dynamical behavior of the ADAM algorithm for nonconvex stochastic optimization[END_REF].

Our approach is inspired by the following continuous-time dynamical system introduced in [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] and referred to as DIN (standing for "dynamical inertial Newton"):

θ(t) Inertial term + α θ(t) Friction term + β ∇ 2 J (θ(t)) θ(t)
Newtonian effects + ∇J (θ(t))

Gravity effect = 0, for t ∈ [0, +∞), ( 1
)
where t is the time parameter which acts as a continuous epoch counter, J is a given loss function (e.g., the empirical loss in DL applications), for now assumed C 2 (twice-differentiable), with gradient ∇J and Hessian ∇ 2 J . It can be shown that solutions to (1) converge to critical Subplots (a-c) represent the trajectories of the parameters θ 1 and θ 2 in R 2 for three choices of hyper-parameters α and β, see (1) for an intuitive explanation. Subplot (d) displays the values of the objective function J (θ 1 , θ 2 ) for the three settings considered.

Neural Networks with Lipschitz Continuous Prediction Function and Loss Functions

We consider DNNs of a general type represented by a function f

: (x, θ) ∈ R M × R P → y ∈ R D
that is locally Lipschitz continuous in θ. This includes for instance feed-forward, convolutional or residual networks used with ReLU, sigmoid, or tanh activation functions. Recall that a function F : R P → R is locally Lipschitz continuous, if for any θ ∈ R P , there exists a neighborhood V of θ and a constant C > 0 such that for any

θ 1 , θ 2 ∈ V , |F (θ 1 ) -F (θ 2 )| ≤ C θ 1 -θ 2 ,
where • is any norm on R P . A function F : R P → R D is locally Lipschitz continuous if each of its coordinates is locally Lipschitz continuous. The variable θ ∈ R P is the parameter of the model (P can be very large), while x ∈ R M and y ∈ R D represent input and output data. For instance, the vector x may embody an image while y is a label explaining its content. Consider further a data set of N samples (x n , y n ) n=1,...,N . Training the network amounts to finding a value of the parameter θ such that, for each input data x n of the data set, the output f (x n , θ) of the model predicts the real value y n with good accuracy. To do so, we follow the traditional approach of minimizing an empirical risk loss function,

R P θ → J (θ) = N n=1 l(f (x n , θ), y n ), (3)
where l : R D × R D → R is a locally Lipschitz continuous dissimilarity measure. In the sequel, for n ≥ 1, we will sometimes denote by J n the n-th term of the sum: J n (θ) l(f (x n , θ), y n ), so that J = N n=1 J n . Despite the non-smoothness and the non-convexity of typical DL loss functions, they generally possess a very strong property sometimes called tameness. We now introduce this notion which is essential to obtain the convergence results of Section 3.

Neural Networks and Tameness in a Nutshell

Tameness refers to a ubiquitous geometrical property of loss functions and constraints encompassing most finite dimensional optimization problems met in practice. Prominent classes of tame objects are piecewise-linear or piecewise-polynomial objects (with finitely many pieces), and more generally, semi-algebraic objects. However, the notion is much more general, as we intend to convey below. The formal definition is given at the end of this subsection (Definition 1).

Informally, sets or functions are called tame when they can be described by a finite number of basic formulas, inequalities, or Boolean operations involving standard functions such as polynomial, exponential, or max functions. We refer to [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF] for illustrations, recipes and examples within a general optimization setting or [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF] for illustrations in the context of neural networks. The reader is referred to [START_REF] Van Den Dries | Tame topology and o-minimal structures[END_REF], [START_REF] Coste | Istituti editoriali e poligrafici internazionali Pisa[END_REF], [START_REF] Shiota | Geometry of subanalytic and semialgebraic sets[END_REF] for foundational material. To apprehend the strength behind tameness it is convenient to remember that it models non-smoothness by confining the study to sets and functions which are union of smooth pieces. This is the so-called stratification property of tame sets and functions. It was this property which motivated the term of tame topology,1 see [START_REF] Van Den Dries | Tame topology and o-minimal structures[END_REF]. In a non-convex optimization settings, the stratification property is crucial to generalize qualitative algorithmic results to non-smooth objects.

Most finite dimensional DL optimization models we are aware of yield tame loss functions J . To understand this assertion and illustrate the wide scope of tameness assumptions, let us provide concrete examples (see also [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF]. Assume that the DNNs under consideration are built from the following traditional components:

-the network architecture describing f is fixed with an arbitrary number of layers of arbitrary dimensions and arbitrary Directed Acyclic Graph (DAG) representing computations, -the activation functions are among classical ones: ReLU, sigmoid, SQNL, RReLU, tanh, APL, soft plus, soft clipping, and many others including multivariate activation functions (norm, sorting), or activation functions defined piecewise with polynomials, exponential and logarithm, -the dissimilarity function l is a standard loss such as p norms, logistic loss or crossentropy, or more generally a function defined piecewise using polynomials, exponentials and logarithms, then one can easily show, by elementary quantifier elimination arguments (property (iii) below), that the corresponding loss, J , is tame.

For the sake of completeness, we provide below the formal definition of tameness and o-minimality.

Definition 1. [o-minimal structure] (Coste, 2000, Definition 1.5) An o-minimal structure on (R, +, .) is a countable collection of sets O = {O q } q≥1 where each O q is itself a collection of subsets of R q , called definable subsets. They must have the following properties, for each q ≥ 1:

(i) (Boolean properties) O q contains the empty set, is stable by finite union, finite intersection and complementation;

(ii) (Lifting property) if A belongs to O q , then A × R and R × A belong to O q+1 .

(iii) (Projection or quantifier elimination property) if Π : R q+1 → R q is the canonical projection onto R q then for any A in O q+1 , the set Π(A) belongs to O q .

(iv) (Semi-algebraicity) O q contains the family of algebraic subsets of R q , that is, every set of the form

{θ ∈ R q | ζ(θ) = 0},
where ζ : R q → R is a polynomial function.

(v) (Minimality property), the elements of O 1 are exactly the finite unions of intervals and points.

A mapping

F : S ⊂ R m → R q is said to be definable in O if its graph is definable in O as a subset of R m × R q .
For illustration of o-minimality in the context of optimization one is referred to [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[END_REF], [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF].

From now on we fix an o-minimal structure O and a set or a mapping definable in O will be called tame.

From DIN to INNA

We describe in this section the construction of our proposed algorithm INNA from the discretization of the second-order ODE (1).

Handling Non-smoothness and Non-convexity

We first show how the formalism offered by Clarke's subdifferential can be applied to generalize (1) to the non-smooth non-convex setting. Recall that the dynamical system (1) is described by,

θ(t) + α θ(t) + β∇ 2 J (θ(t)) θ(t) + ∇J (θ(t)) = 0, (4)
where J is a twice-differentiable potential, α > 0, β > 0 are two hyper-parameters and θ : R + → R P . We cannot exploit (4) directly since in most DL applications J is not twice differentiable (and even not differentiable at all). We first overcome the explicit use of the Hessian matrix ∇ 2 J by introducing an auxiliary variable ψ : R + → R P like in [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF]. Consider the following dynamical system (defined for J merely differentiable),

θ(t) + β∇J (θ(t)) +(α -1 β )θ(t) + 1 β ψ(t) = 0 ψ(t) +(α -1 β )θ(t) + 1 β ψ(t) = 0
, for a.e. t ∈ (0, +∞).

(5)

As explained in [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF], ( 4) is equivalent to (5) when J is twice differentiable. Indeed, one can rewrite (4) into (5) by introducing ψ = -β θβ 2 ∇J (θ) -(αβ -1)θ. Conversely, one can substitute the first line of (5) into the second one to retrieve (4). Note however that (5) does not require the existence of second-order derivatives.

Let us now introduce a new non-convex non-differentiable version of (5). By Rademacher's theorem, locally Lipschitz continuous functions J : R P → R are differentiable almost everywhere. Denote by R the set of points where J is differentiable. Then, R P \ R has zero Lebesgue measure. It follows that for any θ ∈ R P \ R, there exists a sequence of points in R whose limit is this θ . This motivates the introduction of the subdifferential due to [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], defined next.

Definition 2 (Clarke subdifferential of Lipschitz functions). For any locally Lipschitz continuous function F : R P → R, the Clarke subdifferential of F at θ ∈ R P , denoted ∂F (θ), is the set defined by,

∂F (θ) = conv v ∈ R P | ∃(θ k ) k∈N ∈ R N , such that θ k ---→ k→∞ θ and ∇F (θ k ) ---→ k→∞ v , (6)
where conv denotes the convex hull operator. The elements of the Clarke subdifferential are called Clarke subgradients.

The Clarke subdifferential is a nonempty compact convex set. It coincides with the gradient for smooth functions and with the traditional subdifferential for non-smooth convex functions. As already mentioned, and contrarily to the (sub)differential operator, it does not enjoy a sum rule.

Thanks to Definition 2, we can extend (5) to non-differentiable functions. Since ∂J (θ) is a set, we no longer study a differential equation but rather a differential inclusion, given by,

θ(t) + β∂J (θ(t)) +(α -1 β )θ(t) + 1 β ψ(t) 0 ψ(t) +(α -1 β )θ(t) + 1 β ψ(t) 0 , for a.e. t ∈ (0, +∞). (7)
For a given initial condition (θ 0 , ψ 0 ) ∈ R P × R P , we call solution (or trajectory) of this system any absolutely continuous curve (θ, ψ) from R + to R P × R P for which (θ(0), ψ(0)) = (θ 0 , ψ 0 ) and ( 7) holds. We recall that absolute continuity amounts to the fact that θ is differentiable almost everywhere with integrable derivative and,

θ(t) -θ(0) = t 0 θ(s) ds, for t ∈ [0, +∞).
Due to the properties of the Clarke subdifferential, existence of a solution to differential inclusions such as (7) is ensured, see [START_REF] Aubin | Differential inclusions: set-valued maps and viability theory[END_REF]; note however that uniqueness of the solution does not hold in general. We will now use the structure of (7) to build a new algorithm to train DNNs.

Discretization of the Differential Inclusion

To obtain the basic form of our algorithm, we discretize (7) according to the classical explicit Euler method. Given (θ, ψ) a solution of (7) and any time t k , set θ k = θ(t k ) and ψ k = ψ(t k ). Then, at time t k+1 = t k + γ k with γ k positive small, one can approximate θ(t k+1 ) and ψ(t k+1 ) by

θ(t k+1 ) θ k+1 -θ k γ k , ψ(t k+1 ) ψ k+1 -ψ k γ k .
This discretization yields the following algorithm,

       v k ∈ ∂J (θ k ) θ k+1 = θ k + γ k ( 1 β -α)θ k -1 β ψ k -βv k ψ k+1 = ψ k + γ k ( 1 β -α)θ k -1 β ψ k (8)
Although the algorithm above is well-defined for our problem, it is not suited to DL. First, the computation of ∂J (θ k ) is generally not possible since there is no general operational calculus for Clarke's subdifferential; secondly, a mini-batch strategy must be designed to cope with the large dimension of DL problems which makes the absence of sum rule even more critical.

The next section is meant to address these issues and design a practical algorithm.

INNA Algorithm and a New Notion of Steady States

In order to compute or approximate the subdifferential of J at each iteration and to cope with large data sets, J can be approximated by mini-batches, reducing the memory footprint and computational cost of evaluation. For any B ⊂ {1, . . . , N }, let us define

J B : θ → n∈B l(f (x n , θ), y n ). (9)
Unlike in the differentiable case, subgradients do not in general sum up to a subgradient of the sum, that is ∂J B (θ) = n∈B ∂l(f (x n , θ), y n ) in general. To see this, take for example

0 = |•|-|•|, the Clarke subgradient of this function at 0 is {0}, whereas ∂(|0|) + ∂(-|0|) = [-1, 1] + [-1, 1] = [-2, 2].
Standard DL solvers use backpropagation algorithms which implement smooth calculus on non-smooth and non-convex objects. Due to the absence of qualification conditions, the resulting objects are not Clarke subgradients in general. In order to match the real-world practice of DL, we introduce a notion of steady states that corresponds to the stationary points generated by a generic mini-batch approach. As we shall see, this allows both for practical applications and convergence analysis (despite the sum rule failure for Clarke subdifferential). We emphasize once more that our goal is to capture the stationary points that are actually met in practice.

For any B ⊂ {1, . . . , N }, we introduce the following objects,

DJ B = n∈B ∂ [l(f (x n , •), y n )] , DJ = N n=1 ∂ [l(f (x n , •), y n )] . (10)
Observe that, for each B, we have DJ B ⊃ ∂J B and that J B is differentiable almost everywhere with DJ B = ∂J B = {∇J B }, see [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]. In particular DJ = ∂J almost everywhere so that the potential differences with the Clarke subgradient occur on a negligible set. When J is tame the equality even hold on the complement of a finite union of manifolds of dimension strictly lower than P -use the classical stratification results for o-minimal structures, [START_REF] Coste | Istituti editoriali e poligrafici internazionali Pisa[END_REF]. A point satisfying DJ (θ) 0 will be called D-critical. Note that Clarke-critical points (0 ∈ ∂J ) are D-critical points but that the converse is not true. This terminology is motivated by favorable properties: sum and chain rules along curves (see Lemmas 3.2 and 3.3 below) and the existence of a tame Sard's theorem (see Lemma 3.4). To our knowledge, this notion of a steady state has not previously been used in the literature and a direct approach modeling the mini-batch practice has never been considered before.2 While this notion is needed for the theoretical analysis, one should keep in mind that DJ is actually what is computed numerically provided that the automatic differentiation library returns a Clarke subgradient. This computation is usually done with a backpropagation algorithm, similarly to the seminal method of [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF].

Ultimately, one can rewrite (7) by replacing ∂J by DJ , which yields a differential inclusion adapted to study mini-batch approximations of non-smooth loss functions J . This reads,

θ(t) + βDJ (θ(t)) +(α -1 β )θ(t) + 1 β ψ(t) 0 ψ(t) +(α -1 β )θ(t) + 1 β ψ(t) 0 , for a.e. t ∈ (0, +∞). (11)
Discretizing this system gives a workable version of INNA. Let us consider a sequence (B k ) k∈N of nonempty subsets of {1, . . . , N }, chosen independently and uniformly at random with replacement, and a sequence of positive step-sizes (γ k ) k∈N . For a given initialization (θ 0 , ψ 0 ) ∈ R P × R P , at iteration k ≥ 1, our algorithm reads,

(INNA)        v k ∈ DJ B k (θ k ) θ k+1 = θ k + γ k ( 1 β -α)θ k -1 β ψ k -βv k ψ k+1 = ψ k + γ k ( 1 β -α)θ k -1 β ψ k (12)
Here again α > 0 and β > 0 are hyper-parameters of the algorithm. The mini-batch procedure forms a stochastic approximation of the deterministic dynamics obtained by choosing B k = {1, . . . , N }, i.e., when J B k = J (batch version). This can be seen by observing that the vectors v k above may be written as v k = ṽk + η k , where ṽk ∈ DJ (θ k ) and η k compensates for the missing subgradients and can be seen as a zero-mean noise. Hence, INNA admits the following general abstract stochastic formulation,

       w k ∈ DJ (θ k ) θ k+1 = θ k + γ k ( 1 β -α)θ k -1 β ψ k -βw k + ξ k ψ k+1 = ψ k + γ k ( 1 β -α)θ k -1 β ψ k (13)
where (ξ k ) k∈N is a martingale difference noise sequence adapted to the filtration induced by (random) iterates up to k. While ( 12) is the version implemented in practice, its equivalent form ( 13) is more convenient for the convergence analysis of the next section. We stress that the equivalence between ( 12) and ( 13) relies on the use of DJ and would not hold with the use of ∂J as in ( 8).

INNA in its general and practical form is summarized in Table 1.

Inertial Newton Algorithm for Deep Learning (INNA)

Objective function: J = N n=1 J n , with J n : R P → R locally Lipschitz. Hyper-parameters: (α, β) positive. Mini-batches: (B k ) k∈N , nonempty subsets of {1, . . . , N }. step-sizes:

(γ k ) k∈N positive. Initialization: (θ 0 , ψ 0 ) ∈ R P × R P . For k ∈ N:        v k ∈ n∈B k ∂ [J n (θ k )] θ k+1 = θ k + γ k ( 1 β -α)θ k -1 β ψ k -βv k ψ k+1 = ψ k + γ k ( 1 β -α)θ k -1 β ψ k Table 1: INNA in a nutshell.

Convergence Results for INNA

We first state our main result.

Main result: Accumulation Points of INNA are Critical

We now study the convergence of INNA. The main idea here is to prove that the discrete algorithm (12) asymptotically behaves like the solutions of the continuous differential inclusion (11).

In addition to tameness, our main assumption is the following:

Assumption 1 (Stochastic approximation). The sets (B k ) k∈N are taken independently uniformly at random with replacement. The step-size sequence γ k is positive with k γ k = +∞ and satisfies

γ k = o 1 log k , that is lim sup k→+∞ |γ k log k| = 0.
Typical admissible choices are γ k = C(k + 1) -a with a ∈ (0, 1], C > 0. The main theoretical result of the paper follows.

Theorem 3.1 (INNA converges to the set of D-critical points of J ). Assume that for n ∈ {1, . . . , N }, each J n is locally Lipschitz continuous, tame and that the step-sizes satisfy Assumption 1. Set an initial condition (θ 0 , ψ 0 ) and assume that there exists M > 0 such that sup k≥0 (θ k , ψ k ) ≤ M almost surely, where (θ k , ψ k ) k∈N are generated by INNA. Then, almost surely, any accumulation point θ of the sequence (θ k ) k∈N satisfies DJ ( θ) 0. In addition (J (θ k )) k∈N converges.

Before proving Theorem 3.1, we will first make some comments and illustrate this result.

Comments on the Results of Theorem 3.1

-On the step-sizes. First, Assumption 1 offers much more flexibility than the usual O(1/ √ k) assumption commonly used for SGD. We leverage the boundedness assumption on the norms of (θ k , ψ k ), the local Lipschitz continuity and finite-sum structure of J , so that the noise is actually uniformly bounded and hence sub-Gaussian, allowing for much larger step-sizes than in the more common bounded second moment setting. See Benaïm et al. (2005, Remark 1.5) and [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] for more details. The interest of this aggressive strategy is highlighted in Figure 4 of the experimental section.

-On the scope of the theorem. Our result actually holds for general locally Lipschitz continuous tame functions with finite-sum structure and for the general stochastic process under uniformly bounded martingale increment noise. We do not use any other specific structure of DL loss functions. Other variants could be considered depending on the assumptions on the noise, see [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF].

-On D-criticality. The result of Theorem 3.1 states that the bounded discrete trajectories of INNA are attracted by the D-critical points. Recall that D-critical points include local minimizers and thus our theoretical finding agrees with our empirical observations that most initializations lead to "valuable weights" θ and to efficient training. In particular for smooth networks where J is differentiable, limit points of INNA are simply critical points of J . The reader should however remember that when the algorithm is initialized on the D-critical set, the algorithm is stationary as well, even when the initialization is non-Clarke critical. This last point shows that D-points are not introduced to simplify the analysis but to sharply model the use of mini-batch methods on non-convex and non-smooth problems.

Hopefully, in practice one can expect to avoid such points with overwhelming probability. Indeed, [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF] proved that SGD converges with probability one to the set of Clarke-critical points. In other words, D-critical points that are not Clarke-critical are not seen by the dynamics (see also [START_REF] Bianchi | Convergence of constant step stochastic gradient descent for non-smooth non-convex functions[END_REF]. The key argument is to prove that the set of initializations and step-sizes that leads the algorithm to reach points where DJ = ∂J has zero-measure. The same result can be hoped for INNA.

-On the boundedness assumption. The bounded assumption on the iterates is a classical assumption for first or second-order algorithms, see for instance [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF], [START_REF] Duchi | Stochastic methods for composite and weakly convex optimization problems[END_REF]. When using deterministic algorithms (i.e., without mini-batch approximations), properties such as the coercivity of J can be sufficient to remove the boundedness assumption for descent algorithms. This does not remain true when dealing with mini-batch approximations, yet, in the case of INNA, the coercivity of J would guarantee at least that the solutions of the continuous underlying differential inclusion (11) remain bounded. Indeed, we will prove in Section 3.4 that for any solution (θ, ψ) of ( 11), the function

E(θ(t), ψ(t)) 2(1 + αβ)J (θ(t)) + (α -1 β )θ(t) + 1 β ψ(t)
2 is decreasing in time (see Lemma 3.6 hereafter). As a consequence, we cannot have J (θ(t)) ---→ t→∞ ∞ so θ(t) → ∞ due to the coercivity of J . In addition this guarantees ψ(t) → ∞ as well. However, DL loss functions are not coercive in general and studying the boundedness issue in DL or even for non-convex semi-algebraic problems is far beyond the scope of this paper. Let us however mention that it is not uncommon to project the iterates on a given large ball to ensure boundedness; this is a matter for future research.

Preliminary Variational Results

Prior to proving Theorem 3, we extend some results known for the Clarke subdifferential of tame functions to the operator D that we previously introduced. First, we recall a useful result of [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF] which follows from the projection formula in Bolte et al. (2007b).

Lemma 3.2 (Chain rule for the Clarke subdifferential). Let J : R P → R be a locally Lipschitz continuous tame function, then J admits a chain rule, meaning that for all absolutely continuous curves θ : R + → R P , J • θ is differentiable a.e. and for a.e. t ≥ 0,

dJ dt (θ(t)) = θ(t), ∂J (θ(t)) = θ(t), v , ∀v ∈ ∂J (θ(t)). (14)
Note that, even though J is non-differentiable on R P , the function t → J (θ(t)) is differentiable for a.e. t > 0. Indeed, as introduced in Section 2.3.1, an absolutely continuous curve from t ≥ 0 to R P is differentiable for a.e. t > 0. This, combined with the chain-rule of Lemma 3.2 allows differentiating J • θ for a.e. t > 0 whenever J is tame and locally Lipschitz continuous. Additionally, notice that the value of dJ dt (θ(t)) in ( 14) does not depend on the element v taken in ∂J (θ(t)) which justifies the notation θ(t), ∂J (θ(t)) .

Consider now a function with an additive finite-sum structure (such as in DL):

J : R P θ → N n=1 J n (θ), (15)
where each J n : R P → R is locally Lipschitz continuous and tame. We define for any θ ∈ R P

DJ (θ) = N n=1 ∂J n (θ).
The following lemma is a direct generalization of the above chain rule.

Lemma 3.3 (Chain rule for DJ ). Let J be a sum of tame functions like in (15). Let c : [0, 1] → R P be an absolutely continuous curve so that t → J (c(t)) is differentiable almost everywhere. For a.e. t ∈ [0, 1], and for all v ∈ DJ (c(t)),

d dt J (c(t)) = v, ċ(t) .
Proof. By local Lipschitz continuity and absolute continuity, each J n is differentiable almost everywhere and Lemma 3.2 can be applied:

d dt J n (c(t)) = v n , ċ(t) , for all v n ∈ ∂J n (c(t)
) and for a.e. t ≥ 0.

Thus

d dt J (c(t)) = N n=1 d dt J n (c(t)) = N n=1 v n , ċ(t) ,
for any v n ∈ ∂J n (c(t)), for all n = {1, . . . , N }, and for a.e. t ≥ 0. This proves the desired result.

We finish this section with a Sard lemma for D-critical values, in the spirit of Bolte et al. (2007b).

Lemma 3.4 (A Sard's theorem for tame D-critical values). Let,

S = D-crit θ ∈ R P | DJ (θ) 0 , then J (S) is finite.
Proof. The set S is tame and hence it has a finite number of connected components. It is sufficient to prove that J is constant on each connected component of S. Without loss of generality, assume that S is connected and consider θ 0 , θ 1 ∈ S. By Whitney regularity (Van den Dries, 1998, 4.15), there exists a tame continuous path Γ joining θ 0 to θ 1 . Because of the tame nature of the result, we should here conclude with only tame arguments and use the projection formula in Bolte et al. (2007b), but for convenience of readers who are not familiar with this result we use Lemma 3.2. Since Γ is tame, the monotonicity lemma (see for example Kurdyka 1998, Lemma 2) gives the existence of a finite collection of real numbers 0 = a 0 < a 1 < . . . < a q = 1, such that Γ is C 1 on each segment (a j-1 , a j ), j = 1, . . . , q. Applying Lemma 3.2 to each Γ |(a i ,a i+1 ) , we see that J is constant except perhaps on a finite number of points, thus J is constant by continuity.

Proof of Convergence for INNA

Our approach follows the stochastic method for differential inclusions developed in [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF] for which the differential system (11) and its Lyapunov properties play fundamental roles.

The steady states of (11) are given by,

S = (θ, ψ) ∈ R P × R P | 0 ∈ DJ (θ), ψ = (1 -αβ)θ . ( 16
)
These points are initialization values for which the system does not evolve and remains constant.

Observe that the first coordinates of these points are D-critical for J and that conversely any D-critical point of J corresponds to a unique rest point in S.

Definition 3 (Lyapunov function). Let A be a subset of R P × R P , we say that E : R P × R P → R is a Lyapunov function for the set A and the dynamics (11) if, (i) For any solution (θ, ψ) of (11) with initial condition (θ 0 , ψ 0 ) ∈ R P × R P , we have: E(θ(t), ψ(t)) ≤ E(θ 0 , ψ 0 ) a.e. on R.

(ii) For any solution (θ, ψ) of (11) with initial condition (θ 0 , ψ 0 ) ∈ R P × R P \ A, we have: E(θ(t), ψ(t)) < E(θ 0 , ψ 0 ) a.e. on R.

In practice, to establish that a functional is Lyapunov, one can simply use differentiation through chain rule results, with in particular Lemma 3.2. In the context of INNA, we will use Lemma 3.3. To build a Lyapunov function for the dynamics (11) and the set S, consider the two following energy-like functions,

   E min (θ(t), ψ(t)) = (1 - √ αβ) 2 J (θ(t)) + 1 2 (α -1 β )θ(t) + 1 β ψ(t) 2 E max (θ(t), ψ(t)) = (1 + √ αβ) 2 J (θ(t)) + 1 2 (α -1 β )θ(t) + 1 β ψ(t) 2 . ( 17
)
Then the following lemma applies.

Lemma 3.5 (Differentiation along DIN trajectories). Let (θ, ψ) be a solution of (11) with initial condition (θ 0 , ψ 0 ). For a.e. t > 0, θ and ψ are differentiable at t, (11) holds, θ(t)-ψ(t) β ∈ DJ (θ(t)) and

dE min dt (θ(t), ψ(t)) = - √ α θ(t) - 1 √ β ψ(t) -θ(t) 2 dE max dt (θ(t), ψ(t)) = - √ α θ(t) + 1 √ β ψ(t) -θ(t) 2 Proof. Define E λ (θ, ψ) = λJ (θ) + 1 2 (α - 1 β )θ + 1 β ψ 2
. We aim to choose λ so that E λ is a Lyapunov function. Because J is tame and locally Lipschitz continuous, using Lemma 3.3 we know that for any absolutely continuous trajectory θ : R + → R P and for a.e. t > 0,

dJ dt (θ(t)) = θ(t), DJ (θ(t)) = θ(t), v(t) , ∀v(t) ∈ DJ (θ(t)). (18)
Let θ and ψ be solutions of (DIN). For a.e. t ≥ 0, we can differentiate E λ (θ, ψ) to obtain

dE λ dt (θ(t), ψ(t)) =λ θ(t), v(t) + (α - 1 β ) θ(t), (α - 1 β )θ(t) + 1 β ψ(t) + 1 β ψ (t), (α - 1 β )θ(t) + 1 β ψ(t) (19)
for all v(t) ∈ DJ (θ(t)). Using (11), we get 1 β ( θ(t) -ψ(t)) ∈ DJ (θ(t)) and -

ψ(t) = (α -1 β )θ(t) + 1 β ψ(t) a.e. Choosing v(t) = 1 β ( θ(t) -ψ(t)) yields: dE λ dt (θ(t), ψ(t)) = λ θ(t), θ(t) -ψ(t) β -(α - 1 β ) θ(t), ψ(t) - 1 β ψ(t), ψ(t) .
Then, expressing everything as a function of θ and 1 β (ψθ), one can show that a.e. on R + :

dE λ dt (θ, ψ)(t) = -α θ(t) 2 -β θ(t) -ψ(t) β 2 + (λ -αβ -1) θ(t), θ(t) -ψ(t) β = - √ α θ(t) + αβ + 1 -λ 2 √ α θ(t) -ψ(t) β 2 -β - (αβ + 1 -λ) 2 4α θ(t) -ψ(t) β 2 .
We aim to choose λ so that E λ is decreasing that is β -(αβ+1-λ) 2 4α > 0. This holds whenever

λ ∈ (1 - √ αβ) 2 , (1 + √ αβ) 2 . We choose λ min = (1 - √ αβ) 2
, and λ max = (1 + √ αβ) 2 , for these two values we obtain for a.e. t > 0 ,

     Ėλ min (θ(t), ψ(t)) = - √ α θ(t) + 1 √ β θ(t) -ψ(t) 2 Ėλmax (θ(t), ψ(t)) = - √ α θ(t) -1 √ β θ(t) -ψ(t) 2 (20)
Remark finally that by definition E min = E λ min and E max = E λmax .

Recall that S = (θ, ψ) ∈ R P × R P | 0 ∈ DJ(θ), ψ = (1αβ)θ) and define E = E min + E max . By a direct integration argument, we obtain the following lemma.

Lemma 3.6 (E is Lyapunov function for INNA with respect to S). For all (θ 0 , ψ 0 ) / ∈ S and for any solution (θ, ψ) with initial condition (θ 0 , ψ 0 ), E(θ(t), ψ(t)) < E(θ 0 , ψ 0 ), for a.e. t > 0.

(

We are now in position to provide the desired proof.

Proof of Theorem 3.1 Lemmas 3.5 and 3.6 state that E is a Lyapunov function for the set S and the dynamics (11). Let C = {θ ∈ R P | (θ, ψ) ∈ S} which is actually the set of D-critical points of J . Using Lemma 3.4 of Section 3.3, J (C) is finite. Moreover, since E(θ, ψ) = 2(1 + αβ)J (θ) for all (θ, ψ) ∈ S, E takes a finite number of values on S, and in particular, E(S) has empty interior. Denote by L the set of accumulation points of the sequences ((θ k , ψ k )) k∈N produced by (12) starting at (θ 0 , ψ 0 ) and L 1 its projection on R P × {0}. We have the 3 following properties:

-By assumption, we have (θ k , ψ k ) ≤ M almost surely, for all k ∈ N.

-By local Lipschitz continuity ∂J B (θ) is uniformly bounded for θ ≤ M and any B ⊂ {1, . . . , N }, hence the centered noise (ξ k ) k∈N is a uniformly bounded martingale difference sequence.

-By Assumption 1, the sequence (γ k ) k∈N is chosen such that γ k = o( 1 log k ) (see Section 3.2). Then the sufficient conditions of Remark 1.5 of [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF] state that the discrete process (θ k , ψ k ) k∈N asymptotically behaves like the solutions of (11). We can then combine Proposition 3.27 and Theorem 3.6 of [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF], to obtain that the limit set L of the discrete process is contained in the set S where the Lyapunov has vanishing derivatives. Thus, the set L 1 (the set of the first coordinates of all accumulation points) contains only D-critical points of J . In addition, E(L) is a singleton, and for all (θ, ψ) ∈ S, we have E(θ, ψ) = J (θ), so J (L 1 ) is also a singleton and the theorem follows.

Towards Convergence Rates for INNA

In the previous section, connecting INNA to (11) was one of the keys to prove the convergence of the discrete dynamics. Let us now focus on the continuous dynamical system (7) in the deterministic case where J and ∂J are not approximated anymore-we thus no longer use DJ although this would be possible but would require more technical proofs. In this section and in this section only, we pertain to loss functions J that are real semi-algebraic (a particular case of tame functions).3 Recall that a set is called semi-algebraic if it is a finite union of sets of the form,

{θ ∈ R P | ζ(θ) = 0, ζ i (θ) < 0}
where ζ, ζ i are real polynomial functions. A function is called semi-algebraic if its graph is semi-algebraic.

We will characterize the convergence rate of the solutions of the continuous-time system (7) to critical points. Let us first introduce an essential mechanism to obtain such convergence rates: the Kurdyka-Łojasiewicz (KL) property.

The Non-smooth Kurdyka-Łojasiewicz Property for the Clarke Subdifferential

The non-smooth Kurdyka-Łojasiewicz (KL) property, as introduced in [START_REF] Bolte | Characterizations of Łojasiewicz inequalities: subgradient flows, talweg, convexity[END_REF], is a measure of "amenability to sharpness" (as illustrated at the end of Section 4.3). Here we provide a uniform version for the Clarke subdifferential of semi-algebraic functions as in Bolte et al. (2007b) and [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. In the sequel we denote by "dist" any given distance on R P .

Lemma 4.1 (Uniform Non-smooth KL Property for the Clarke Subdifferential). Let K be a nonempty compact set and let L : R P → R be a semi-algebraic locally Lipschitz continuous function. Assume that L is constant on K, with value L . Then there exist ε > 0, δ > 0, a ∈ (0, 1) and ρ > 0 such that, for all

v ∈ v ∈ R P | dist(v, K) < ε ∩ v ∈ R P | L < L(v) < L + δ , it holds that, ρ(1 -a) (L(v) -L(v)) -a dist (0, ∂L(v)) > 1. ( 22
)
The proof directly follows from the general inequality provided in Bolte et al. (2007b) or the local result of Bolte et al. (2007b) with the compactness arguments of Bolte et al. (2014, Lemma 6). In the sequel, we make an abuse of notation by writing ∂J (•) dist(0, ∂J (•)). To obtain a convergence rate we will use inequality ( 22) on the Lyapunov function E. But first we state a general result of convergence that is built around the KL property.

A General Asymptotic Rate Result

We state a general theorem that leads to the existence of a convergence rate. This theorem will hold in particular for (7). We start by stating the result.

Theorem 4.2. Let X : [0, +∞) → R P be a bounded absolutely continuous trajectory and let L : R P → R be a semi-algebraic locally Lipschitz continuous function. If there exists c 1 > 0 such that for a.e. t > 0, dL dt (X(t)) ≤ -c 1 (∂L)(X(t)) 2 , (i)

then L(X(t)) converges to a limit value L and,

|L(X(t)) -L | = O 1 t .
If in addition there exists c 2 > 0 such that for a.e. t > 0,

c 2 Ẋ(t) ≤ (∂L)(X(t)) , (ii)
then, X converges to a critical point of L with a rate of the form O(1/t b ) with b > 0.4

Proof. We first prove the convergence of L(X(•)). Suppose that (i) holds. Since X is bounded and L is continuous, L(X(•)) is bounded. Moreover, from (i), L(X(•)) is decreasing, so it converges to some value L . To simplify suppose L ≥ 0 and L = 0. Define,

I = x ∈ R P | L(x) = 0 .
Suppose first that there exists s ≥ 0, such that X(s) ∈ I. Since L(X(•)) is decreasing with limit 0, then for all t ≥ s, L(X(t)) = 0 and the convergence rate holds true. Let us thus assume that for all t ≥ 0, L(X(t)) > 0. The trajectory X is bounded in R P , hence there exists a compact set C ⊂ R P such that X(t) ∈ C for all t ≥ 0. Define K = I ∩ C. It is a compact set since I is closed (by continuity of L) and C is compact. Moreover, L is constant on K. As such by Lemma 4.1, there exist ε > 0, δ > 0, a ∈ (0, 1) and a constant ρ > 0 such that for all

v ∈ v ∈ R P , dist(v, K) < ε ∩ {0 < L(v) < δ} , it holds that ρ(1 -a) (L(v)) -a dist (0, ∂L(v)) > 1.
We have L(X(t)) → 0 so there exists t 0 ≥ 0 such that for all t ≥ t 0 , 0 < L(X(t)) < δ. Without loss of generality, we assume t 0 = 0. Similarly, we have dist(X(t), K) → 0, so we may assume that for all t ≥ 0, dist (X(t), K) < ε. Thus, for all t ≥ 0,

ρ(1 -a)L(X(t)) -a ∂L(X(t)) > 1.
Going back to assumption (i), for a.e. t > 0, one has

dL dt (X(t)) ≤ -c 1 (∂L)(X(t)) 2 ,
but the KL property implies that for a.e. t > 0,

-∂L(X(t)) 2 < - 1 ρ 2 (1 -a) 2 L(X(t)) 2a . Therefore, dL dt (X(t)) < - c 1 ρ 2 (1 -a) 2 L(X(t)) 2a .
We consider two cases depending on the value of a. If 0 < a ≤ 1/2, then for t large, L(X(t)) < 1 so -L(X(t)) 2a < -L(X(t)) and hence,

dL dt (X(t)) < - c 1 ρ 2 (1 -a) 2 L(X(t)),
so we obtain a linear rate. When 1/2 < a < 1, we have for a.e. t > 0,

L(X(t)) -2a d dt L(X(t)) = 1 1 -2a d dt L(X(t)) 1-2a < - c 1 (ρ 2 (1 -a) 2 ) , (23)
with 1 -2a < 0. We can integrate from 0 to t > 0:

L(X(t)) 1-2a > (2a -1)c 1 ρ 2 (1 -a) 2 t + L(X(0)) 1-2a > (2a -1)c 1 ρ 2 (1 -a) 2 t.
Since 1 1-2a < -1, one obtains a convergence rate of the form O t 1 1-2a . In both cases the rate is at least O 1 t . We assume now that both (i) and (ii) holds and prove convergence of the trajectory with a convergence rate. Let t > s > 0, by the fundamental theorem of calculus and the triangular inequality,

X(t) -X(s) ≤ t s Ẋ(τ ) dτ ≤ t s Ẋ(τ ) dτ. ( 24
)
We wish to bound Ẋ using L. Using the chain rule (Lemma 3.2 of Section 3.3), for a.e. τ > 0,

d dτ L(X(τ )) 1-a = (1 -a)L(X(τ )) -a Ẋ(τ ), (∂L)(X(τ )) . (25)
Then, from (i), we deduce that for a.e. τ > 0,

Ẋ(τ ), (∂L)(X(τ )) = dL dτ (X(τ )) ≤ -c 1 (∂L)(X(τ )) 2 , ( 26
) so d dτ L(X(τ )) 1-a ≤ -c 1 (1 -a)L(X(τ )) -a (∂L)(X(τ )) 2 . ( 27
)
The KL property ( 22) implies that for a.e. τ > 0,

-(1 -a)L(X(τ )) -a (∂L)(X(τ )) < - 1 ρ . ( 28
)
Putting this in ( 27) and using assumption (ii) we finally obtain

d dt L(X(τ )) 1-a < - c 1 ρ (∂L)(X(τ )) ≤ - c 1 c 2 ρ Ẋ(τ ) . ( 29
)
We can use that in (24),

X(t) -X(s) ≤ - ρ c 1 c 2 t s d dt L(X(τ )) 1-a dτ = ρ c 1 c 2 (L(X(s)) 1-a -L(X(t)) 1-a ).
(30)

Then, using the convergence rate we already proved for L, we deduce that the Cauchy criterion holds for X inside the compact (hence complete) subset C ⊂ R P containing the trajectory. Thus, X converges, and from (i) we have that lim inf t→+∞ ∂L(X(t)) = 0 because ∂L has closed graph. This shows that the limit is a critical point of L. Finally, taking the limit in (30) and using the convergence rate of L we obtain a rate for X as well.

Remark 1. Theorem 4.2 takes the form of a general recipe to obtain a convergence rate since it may be applied in many cases, to curves or flows, provided that a convenient Lyapunov function is given. Note also that it is sufficient for assumptions (i) and (ii) to hold only after some time t 0 > 0 as in such case, one could simply do a time shift to use the theorem.

Application to INNA

We now apply Theorem 4.2 to the deterministic continuous dynamical model of INNA (7).

Theorem 4.3 (Convergence rates). Suppose that J is semi-algebraic locally Lipschitz continuous and lower bounded. Then, any bounded trajectory (θ, ψ) that solves (7) converges to a point ( θ, ψ) ∈ S, with a convergence rate of the form O t -b with b > 0. Moreover, J (θ(t)) converges to its limit J with rate J (θ(t)) -J = O 1 t . Proof. Let (θ, ψ) be a bounded solution of (7). We would like to use Theorem 4.2 with X = (θ, ψ), and a well-chosen function. In the proof of Theorem 3.1 we proved a descent property along the trajectory for the function

E(θ, ψ) = 2(1 + αβ)J (θ) + (α -1 β )θ + 1 β ψ 2
. This function is semi-algebraic, locally Lipschitz continuous, so it remains to prove that (i) and (ii) hold for E along (θ, ψ).

For t ≥ 0, denote w(t) = (α -1 β )θ(t) + 1 β ψ(t), then according to Lemma 3.5 for a.e. t > 0,

dE dt (θ(t), ψ(t)) = - √ α θ(t) - 1 √ β ψ(t) -θ(t) 2 - √ α θ(t) + 1 √ β ψ(t) -θ(t) 2 = -2α θ(t) 2 - 2 β ψ(t) -θ(t) 2 = -2α θ(t) 2 - 2 β β∂J (θ(t)) 2 = -2α -β∂J (θ(t)) -w(t) 2 -2β ∂J (θ(t)) 2 . (31)
On the other hand, by standard results on the sum rule, we have for all (θ, ψ) ∈ R P × R P ,

∂E(θ, ψ) = 2   (1 + αβ)∂J (θ) + (α -1 β ) (α -1 β )θ + 1 β ψ 1 β (α -1 β )θ + 1 β ψ   , (32)
so for a.e. t > 0,

∂E(θ(t), ψ(t)) 2 4 = (1 + αβ)∂J (θ(t)) + (α - 1 β )w(t) 2 + 1 β w(t) 2 . ( 33
)
We wish to find c 1 > 0, such that 1 2 dE dt + c 1 4 ∂E 2 < 0. This follows from the following claim. Claim: let r 1 > 0, r 2 ∈ R, r 3 > 0, then there exist C 1 and C 2 two positive constants such that for any a, b ∈ R,

C 1 (a 2 + b 2 ) ≤ (r 1 a + r 2 b) 2 + r 3 b 2 ≤ C 2 (a 2 + b 2 ). (34)
Indeed, the function Q : (a, b) → (r 1 a + r 2 b) 2 + r 3 b 2 is a positive definite quadratic form, C 1 and C 2 can be taken to be two eigenvalues of the positive definite matrix which represents Q. Hence (34) holds for all a and b.

Applying the previous claim to ( 33) and ( 31) leads to the existence of c 1 > 0 such that for a.e.

t > 0, dE dt (θ(t), ψ(t)) ≤ -c 1 ∂E(θ(t), ψ(t)) 2 ,
so assumption (i) holds for INNA. It now remains to show that (ii) of Theorem 4.2 holds i.e., that there exists c 2 > 0 such that for (θ, ψ) solution of ( 7) and for a.e. t > 0, ∂E(θ(t), ψ(t)) 2 ≥ c 2 θ(t) 2 + ψ(t) 2 . Using ( 7) and ( 33) we obtain:

∂E(θ(t), ψ(t)) 2 4 = 1 β (1 + αβ) θ(t) + (α - 1 β ) - 1 β (1 + αβ) ψ(t) 2 + 1 β 2 ψ(t) 2 , ( 35
)
and applying the claim (34) again to (35) one can show that there exist c 2 > 0, such that for a.e. t > 0,

∂E(θ(t), ψ(t)) 2 ≥ c 2 θ(t) 2 + ψ(t) 2 .
So assumption (ii) holds for (7). To conclude, we can apply Theorem 4.2 to (7) and the proof is complete.

Remark 2. (a) Since the discrete algorithm INNA asymptotically resembles its continuous-time version (see the proof of Theorem 3.1), the results above suggest that similar behaviors and rates could be hoped for INNA itself. Yet, these results remain difficult to obtain in the case of DL, in particular in the mini-batch setting because of the noise (ξ k ) k∈N . (b) The proof above is significantly simpler when αβ > 1 since [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] proved that in this case, ( 7) is equivalent to a gradient system, thus assumptions (i) and (ii) of Theorem 4.2 instantly hold. (c) Theorems 4.2 and 4.3 can be adapted to the case where the Clarke subdifferential is replaced by DJ , but we do not state it here for the sake of simplicity. (d) Theorems 4.2 and 4.3 are actually valid by assuming that J belongs to a polynomially bounded o-minimal structure. One of the most common instance of such structures is the one given by globally subanalytic sets (as illustrated in a example below). We refer to Bolte et al. (2007a) for a definition and further references.

Let us now comment the results of Theorem 4.3. First, we restrained here to semi-algebraic loss functions J , which are a subset of tame loss functions. Most networks, activation functions and dissimilarity measures mentioned in Section 2.2 fall into this category. Nonetheless, the loss functions of the DL experiments of Section 5.2 are not semi-algebraic. Indeed, the dissimilarity measure l used is the cross-entropy:

l(f (x n , θ), y n ) = -D d=1 1 [yn] d =1 log([f (x n , θ)] d ).
Such a function cannot be described by polynomials and presents a singularity whenever [f (x n , θ)] d = 0. Fortunately, due to the numerical precision but also to the "soft-max" functions often used in classification experiments, the outputs of the network f , for inputs restricted to a compact set, have values in [ε, 1] for some small ε > 0. Therefore, the singularity at 0 is harmless and the cross-entropy acts as a globally subanalytic function. As a consequence the non-smooth Łojasiewicz inequality holds, and the theorems apply (see also numerical experiments).

The rate of convergence of the trajectory in Theorem 4.3 is non-explicit in the sense that the exponent b > 0 is unknown in general. In the light of the proof of Theorem 4.2, this exponent depends on the KL exponent a of the Lyapunov function, which is itself hard to determine in practice. However, the intuition is that small exponents a may yield faster convergence rates (indeed, when a ∈ (0, 1/2) we actually have a linear rate). As an example, for the function: t ∈ R → |t| c with c > 1, the exponent at 0 is a = 1 -1 c and thus, the closer c is to 1, the smaller a is, and the faster the convergence becomes.

Experiments

In this section we first discuss the role and influence of the hyper-parameters of INNA using the 2D example given in Figure 1. We then compare INNA with SGD, ADAGRAD and ADAM on deep learning problems for image recognition.

Understanding the Role of the Hyper-parameters of INNA

Both hyper-parameters α and β can be seen as damping coefficients from the viewpoint of mechanics as discussed by [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] and sketched in the introduction. Recall the second-order time continuous dynamics which served as a model to the design of INNA:

θ(t) + α θ(t) + β ∇ 2 J (θ(t)) θ(t) + ∇J (θ(t)) = 0.
This differential equation was inspired by Newton's second law of dynamics asserting that the acceleration of a material point coincides with the sum of forces applied to the particle. As recalled in the introduction three forces are at stake: the gravity and two friction terms. The parameter α calibrates the viscous damping intensity as in the Heavy Ball friction method of [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF]. It acts as a dissipation term but it can also be seen as a proximity parameter of the system with the usual gradient descent: the higher α is, the more DIN behaves like a pure gradient descent. 5 On the other hand the parameter β can be seen as a Newton damping which takes into account the geometry of the landscape to brake or accelerate the dynamics in an adaptive anisotropic fashion, see [START_REF] Alvarez | A dynamical system associated with Newton's method for parametric approximations of convex minimization problems[END_REF], [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] for further insights.

We now turn our attention to INNA, and illustrate the versatility of the hyper-parameters α and β in this case. We proceed on a 2D visual non-smooth ill-conditioned example à la Rosenbrock, see Figure 1. For this example, we aim to find the minimum of the function J (θ 1 , θ 2 ) = 100(θ 2 -|θ 1 |) 2 + |1θ 1 |. This function has a V-shaped valley, and a unique critical point at (1, 1) which is also the global minimum. Starting from the point (-1, 1.5) (the black cross), we apply INNA with constant steps γ k = 10 -4 . Figure 1 shows that when β is too small, the trajectory presents many transverse oscillations as well as longitudinal ones close to the critical point (subplot a). Then, increasing β significantly reduces transverse oscillations (subplot b). Finally, the longitudinal oscillations are reduced by choosing a higher α (subplot c). In addition, these behaviors are also reflected in the values of the objective function (subplot d). The orange curve (first setting) presents large oscillations. Moreover, looking at the red curve, corresponding to plot (c), there is a short period between 20, 000 and 60, 000 iterations when the decrease is slower than for the other values of α and β, but still it presents fewer oscillations. In the longer term, the third choice (α = 1.3, β = 0.1) provides remarkably good performance.

The choice of these hyper-parameters may come with rates of convergence for convex and strongly convex smooth functions [START_REF] Attouch | First-order optimization algorithms via inertial systems with Hessian driven damping[END_REF]. Following this work, one may also consider to make α and β vary in time (for example like the famous Nesterov damping coefficient α t ). In our DL experiments we will however keep these parameters constant so that our theorems still hold. Yet, different behaviors depending on (α, β) can also be observed for DL problems as illustrated on Figure 2 and described next. Although we did not evidence some universal method to choose (α, β), we used mechanical intuitions to tune these parameters. The coefficient α induces viscous damping, thus one may try to reduce it when convergence appears to be slow. On the other hand, one may want to increase β when large oscillations are observed. Yet, since β affects directly the subgradient effect in (12), taking β too large may jeopardize the numerical stability of the algorithm. It would be interesting to relate the roles of these coefficients with more structural aspects of the dynamics. Indeed, as mentioned in Remark 2-b, when αβ ≥ 1, (11) can be shown to be a gradient system. On the other hand, when αβ < 1 the dynamics is of a different type, we refer to [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping: Application to optimization and mechanics[END_REF] for further comments on this matter.

Training a DNN with INNA

Before comparing INNA to state-of-the-art algorithms in DL, we first describe the methodology that we followed.

Methodology

-We train a DNN for classification using the three most common image data sets (MNIST, CIFAR-10, CIFAR-100) [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF]. These data sets are composed of 60, 000 small images associated with a label (numbers, objects, animals, etc.). We split the data sets into 50, 000 images for training and 10, 000 for testing.

-Regarding the network, we use a slightly modified version of the popular Network in Network (NiN) [START_REF] Lin | Network in Network[END_REF]. It is a reasonably large convolutional network with P ∼ 10 6 parameters to optimize. We use ReLU activation functions.

-The dissimilarity measure l that is used in the empirical loss J given by ( 3) is set to the cross-entropy. The loss function J is optimized with respect to θ (the weights of the DNN) on the training data. The classification accuracy of the trained DNN is measured using the test data of 10, 000 images. Measuring the accuracy boils down to counting how many of the 10, 000 were correctly classified (in percentage).

-Based on the results of Section 5.1, we run INNA for four different values of (α, β):

(α, β) ∈ {(0.1, 0.1), (0.5, 0.1), (0.5, 0.5), (0.5, 1)} .

Given an initialization of the weights θ 0 , we initialize ψ 0 such that the initial velocity is in the direction of -∇J (θ 0 ). More precisely, we use

ψ 0 = (1 -αβ)θ 0 -(β 2 -β)∇J (θ 0 ).
-We compare our algorithm INNA with the classical SGD algorithm and the popular ADAGRAD [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF] and ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] algorithms. At each iteration k, we compute the approximation of ∂J (θ) on a subset B k ⊂ {1, . . . , 50, 000} of size 32. The algorithms are initialized with the same random weights (drawn from a normal distribution). Five random initializations are considered for each experiment.

-Regarding the selection of step-sizes, ADAGRAD and ADAM both use an adaptive procedure based on past gradients, see [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF], [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. For the other two algorithms (INNA and SGD), we use the classical step-size schedule γ k = γ 0 √ k+1 , which meets Assumption 1. For all four algorithms, choosing the right initial step length γ 0 is often critical in terms of efficiency. We choose this γ 0 using a grid-search: for each algorithm we select the initial step-size that most decreases the training error J after fifteen epochs (one epoch consisting in a complete pass over the data). The test data is not used to choose the initial step-size nor other hyper-parameters. Note that we could use more flexible step-size schedules but chose a standard schedule for simplicity. Other decay schemes are considered in Figure 4.

For these experiments, we used Keras 2.2.4 [START_REF] Chollet | Keras[END_REF] with Tensorflow 1.13.1 [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF] as backend. The INNA algorithm is available in Pytorch, Keras and Tensorflow: https://github.com/camcastera/Inna-for-DeepLearning/ [START_REF] Castera | INNA for deep learning[END_REF].

Results

Figure 2 displays the training loss J and test accuracy with respect to epochs for INNA in its four hyper-parameter configurations considered and for the three data sets considered. Figure 3 displays the performance of INNA with the hyper-parameter configuration that led to the smallest average training error in Figure 2, with comparison to SGD, ADAGRAD and ADAM. In these two figures (and also in subsequent Figure 4), solid lines represent mean values and pale surfaces represent the best and worst runs in terms of training loss and validation accuracy over five random initializations.

Figure 2 suggests that the tuning of the hyper-parameters α and β is not crucial to obtain satisfactory results both for training and testing. It mostly affects the training speed. Thus, INNA looks quite stable with respect to these hyper-parameters. Setting (α, β) = (0.5, 0.1) appears to be a good default choice when necessary. Nevertheless, tuning these hyper-parameters is of course advised to get the most from INNA.

Figure 3 shows that best performing methods achieve state-of-the art accuracy using NiN and represent what can be achieved with a moderately large network and coarse grid-search Finally, let us point out that although ADAM was faster in the experiments of Figure 3, INNA can outperform ADAM using the slow step-size decay discussed in Section 3.2. Indeed, in the previous experiments we used a standard decreasing step-size of the form γ 0 / √ k + 1 for simplicity, but Assumption 1 allows for step-sizes decreasing much more slowly. As such, we also considered decays of the form γ 0 (k + 1) -q with q ≤ 1/2. The results are displayed on top of Figure 4. Except when q is too small (too slow decay, e.g., q = 1/16), these results show that some decays slower than q = 1/2 make INNA a little faster than any of the other algorithms we tried. In particular, with a step-size decay proportional to k -1/4 , INNA outperforms ADAM (bottom of Figure 4). This suggests that tuning q can also significantly accelerate the training process.

Conclusion

We introduced a novel stochastic optimization algorithm featuring inertial and Newtonian behavior motivated by applications to deep learning. We provided a powerful algorithmic convergence In the legend, k -q means a step-size decay at iteration k of the form γ k = γ 0 k -q . The bottom row show the comparison between INNA with a well-chosen step-size decay and ADAM.

analysis under weak hypotheses applicable to most DL problems. We also provided new general results to study differential inclusions on Clarke subdifferential and obtain convergence rates for the continuous-time counterpart of our algorithm. We would like to point out that, apart from SGD [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF], the convergence of concurrent methods in such a general setting is still an open question. Our result seems moreover to be the first one to be able to rigorously handle the analysis of mini-batch sub-sampling for ReLU DNNs via the introduction of the D-critical points. Our experiments show that INNA is very competitive with state-of-the-art algorithms for DL but also very malleable. We stress that these numerical manipulations were performed on substantial DL benchmarks with only limited algorithm tuning. This facilitates reproducibility and allows staying as close as possible to the reality of DL applications in machine learning.

Figure 1 :

 1 Figure 1: Illustration of INNA applied to the non-smooth function J(θ 1 , θ 2 ) = 100(θ 2 -|θ 1 |) 2 + |1θ 1 |. Subplots (a-c) represent the trajectories of the parameters θ 1 and θ 2 in R 2 for three choices of hyper-parameters α and β, see (1) for an intuitive explanation. Subplot (d) displays the values of the objective function J (θ 1 , θ 2 ) for the three settings considered.

Figure 2 :

 2 Figure 2: Analysis of the sensibility of INNA to the choice of α and β using NiN for three different image classification problems. Top: logarithm of the loss function J (θ) during the training. Bottom: classification accuracy on the test set.

Figure 3 :

 3 Figure 3: Comparison of INNA with state-of-the-art algorithms SGD, ADAM and ADAGRAD. Top: logarithm of the loss function J (θ) during the training. Bottom: classification accuracy on the test set.

Figure 4 :

 4 Figure 4: On top: Training loss of INNA on three image classification problems with various step-size decays.In the legend, k -q means a step-size decay at iteration k of the form γ k = γ 0 k -q . The bottom row show the comparison between INNA with a well-chosen step-size decay and ADAM.

"La topologie modérée" wished for by Grothendieck.

Following the first arXiv preprint of this work[START_REF] Castera | An inertial Newton algorithm for deep learning[END_REF],Bolte and Pauwels (2020a) have further developed the present ideas and in particular the connection to the backpropagation algorithm.

We could extend the results of this section to more general objects including analytic functions on bounded sets. The semi-algebraicity assumption is made here for the sake of clarity.

In some cases we even have linear rates or finite convergence as detailed in the proof.

This is easier to see when one rescales J by α.

Acknowledgments

The authors acknowledge the support of the European Research Council (ERC FACTORY-CoG-6681839), the Agence Nationale de la Recherche (ANR 3IA-ANITI, ANR-17-EURE-0010 CHESS, ANR-19-CE23-0017 MASDOL) and the Air Force Office of Scientific Research (FA9550-18-1-0226).

Part of the numerical experiments were done using the OSIRIM platform of IRIT, supported by the CNRS, the FEDER, Région Occitanie and the French government (http://osirim. irit.fr/site/en). We thank the development teams of the following libraries that were used in the experiments: Python [START_REF] Rossum | Python reference manual[END_REF], Numpy [START_REF] Walt | The NumPy array: a structure for efficient numerical computation[END_REF], Matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF], Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF], Tensorflow and Keras [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF][START_REF] Chollet | Keras[END_REF].

The authors thank the anonymous reviewers for their comments which helped to improve the paper and thank Hedy Attouch and Sixin Zhang for useful discussions.