Jérôme Bolte

Edouard Pauwels

A mathematical model for automatic differentiation in machine learning

Automatic differentiation, as implemented today, does not have a simple mathematical model adapted to the needs of modern machine learning. In this work we articulate the relationships between differentiation of programs as implemented in practice and differentiation of nonsmooth functions. To this end we provide a simple class of functions, a nonsmooth calculus, and show how they apply to stochastic approximation methods. We also evidence the issue of artificial critical points created by algorithmic differentiation and show how usual methods avoid these points with probability one.

Introduction

Optimization algorithms based on backpropagation oracles, and more generally on automatic or algorithmic differentiation (AD) [START_REF] Speelpenning | Compiling fast partial derivatives of functions given by algorithms[END_REF][START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF], are one of the most widely used training tools for modern learning architectures [START_REF] Bottou | The tradeoffs of large scale learning[END_REF][START_REF] Lecun | Deep learning[END_REF][START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Chizat | On the global convergence of gradient descent for overparameterized models using optimal transport[END_REF][START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF][START_REF] Baydin | Automatic differentiation in machine learning: a survey[END_REF][START_REF] Castera | An Inertial Newton Algorithm for Deep Learning[END_REF]. They often rely on popular numerical implementations as TensorFlow or PyTorch [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF][START_REF] Paszke | Automatic differentiation in pytorch[END_REF]. However, for nonsmooth, nonconvex losses, AD does not have a stable theory [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF][START_REF] Griewank | On Lipschitz optimization based on gray-box piecewise linearization[END_REF][START_REF] Barton | Computationally relevant generalized derivatives: theory, evaluation and applications[END_REF][START_REF] Kakade | Provably correct automatic sub-differentiation for qualified programs[END_REF][START_REF] Griewank | Treating artificial neural net training as a nonsmooth global optimization problem[END_REF][START_REF] Griewank | Beyond the Oracle: Opportunities of Piecewise Differentiation[END_REF][START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF], matching the actual practice. We wish to present a simple mathematical framework addressing this issue. Let us progressively explain our approach.

What is backpropagation?

Algorithmic differentiation acts on programs not on functions: To convey this fact we carry out a small experiment in TensorFlow [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF] with the function relu : t Þ Ñ maxt0, tu, see Appendix A.2 for implementation details. Algorithmic differentiation is displayed in Figure 1, in particular, we have relu 1 p0q " 0. Consider the two functions

relu 2 : t Þ Ñ relup´tq `t, relu 3 : t Þ Ñ 1 2 preluptq `relu 2 ptqq.
As mathematical functions on R these are equal to relu. However TensorFlow returns relu 1 2 p0q " 1 and relu 1 3 p0q " 1{2 (Figure 1). Indeed, AD does not act on functions, but on their representations, i.e., on programs. Different programs implementing the same function may provide different results, beyond numerical precision; we refer to this as the spurious behaviour of AD for nonsmooth functions 2 . Let us explore this phenomenon further. The function zero : t Þ Ñ relu 2 ptq ´reluptq, outputs constantly 0 but AD gives zero 1 p0q " 1. More generally, one can modify the value of the derivative of a given function at prescribed arguments (Figure 1). This may generate artificial critical points; for instance x Ñ x ´zero is the identity but its derivative at 0 according to AD is 0.

˚Authors in alphabetical order. 2 The validity domain of AD is restricted in theory to smooth functions [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF], yet it is common practice to use it for nonsmooth functions. This discussion was limited to univariate functions, but these pathologies grow in size and in complexity when occurring in higher dimensions. Besides, as the "compositional depth" of functions increases the phenomenon gets more complex, making the geometry of artificial point difficult to grasp.

Canonical surjection between programs functions: Numerical programs combine basic mathematical functions within an algorithm and return an output. This can be understood in two ways:

• Computer science: it is a sequence of instructions with numerical inputs-outputs,

• Mathematics: the program is a function3 of its arguments.

It is tempting to identify both, but functions can be represented by different programs. This defines a surjection F mapping a program to a function (in the class of functions "accessible through coding").

Algorithmic differentiation: As presented above, AD is an operation on programs, A which takes as argument a program and returns a program with the same input variables. This operation can be "pushed" to the space of functions using the canonical surjection F. Remarkably, if we restrict ourselves to programs P which only smoothly combine smooth functions, then we have the following fundamental relation, depicted in Figure 2:

FpApPqq " ∇FpPq. (1)
In other words, algorithmic differentiation of a program which smoothly combines smooth functions, is equivalent, through the canonical surjection, to derivation. However practitioners use AD and backpropagation beyond smooth programs with nonsmooth elementary functions or program branching for instance. Can we find a proper operational interpretation of this widespread practice?

Algorithmic differentiation cannot be represented through a variational operator At first, it is tempting to simply use AD to induce a differential operator on functions generalizing classical differentiation. This operator, say B A , should:

(a) encompass the outputs of algorithmic differentation for all functions (b) be such that 0 is an element of B A preluq at 0.

Unfortunately such an operator does not exist:

Theorem 1 (Algorithmic differentiation does not induce an operator on functions) There is no nontrivial operator on functions satisfying paq and pbq.

Contribution and related work

We address this impossibility result and provide a class of functions together with an operational nonsmooth differential calculus which is able to cope with spurious behaviours.

Elementary selections and selection derivatives:

We introduce a new class of nonsmooth nonconvex functions, encompassing most objective functions met in machine learning, having appealing stability properties. This allows us to define simple differential objects called selection derivatives. Selection derivatives turn out to have an operational calculus adapted to the analysis of many learning methods, as backpropagation or stochastic first order methods. They thus provide an operational model to capture nonsmooth AD as implemented in current numerical software.

Algorithmic differentiation, algorithms

This framework allows to formalize properly the relationships between, functions, algorithmic differentiation and capture the corresponding notion of critical points as met in practice. These characterize the set of attractors (limit points) for stochastic approximation algorithms based on nonsmooth backpropagation [START_REF] Robbins | A stochastic approximation method[END_REF][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]. It is important to stress that these attractors, which models sharply the whole scope of AD-induced stationnarity, are different from the traditional notions as Clarke criticality [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF][START_REF] Rockafellar | Variational analysis[END_REF][START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF]. This is described in Theorems 3 and 4.

Avoidance of traps: As sketched above and in the introduction AD produces artificial critical points, i.e. stationary points which are not Clarke critical. These points have a parasitic nature which could be detrimental to training purposes, were they met in practice. We show that randomly initialized mini-batch stochastic gradient method do not lead to artificial critical points (Theorem 4). This result applies to modern machine learning software libraries based on AD [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF][START_REF] Paszke | Automatic differentiation in pytorch[END_REF], seen as performing operation over the reals, without any modification. Although AD may have unpredictable behavior in nonsmooth contexts, both theoretically and numerically, this result justifies theoretically that the practical impact is somewhat negligible in the context of common machine learning usage.

Related work: Spurious behaviour of AD in nonsmooth context has been investigated in [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF][START_REF] Griewank | On Lipschitz optimization based on gray-box piecewise linearization[END_REF][START_REF] Griewank | First-and second-order optimality conditions for piecewise smooth objective functions[END_REF][START_REF] Barton | Computationally relevant generalized derivatives: theory, evaluation and applications[END_REF][START_REF] Kakade | Provably correct automatic sub-differentiation for qualified programs[END_REF][START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF]. In particular, [START_REF] Griewank | First-and second-order optimality conditions for piecewise smooth objective functions[END_REF][START_REF] Kakade | Provably correct automatic sub-differentiation for qualified programs[END_REF] considers qualification conditions allowing to construct AD algorithms which compute proper Clarke subgradients [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF][START_REF] Rockafellar | Variational analysis[END_REF][START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF]. However qualification is extremely hard to check and almost impossible to enforce in practice. Let us also mention [START_REF] Barton | Computationally relevant generalized derivatives: theory, evaluation and applications[END_REF] which uses the notion of lexicographic derivatives, but, at this day, algorithmic computations are limited to forward mode for the moment which is of little use in machine learning. [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF][START_REF] Griewank | First-and second-order optimality conditions for piecewise smooth objective functions[END_REF][START_REF] Griewank | On Lipschitz optimization based on gray-box piecewise linearization[END_REF][START_REF] Griewank | Treating artificial neural net training as a nonsmooth global optimization problem[END_REF][START_REF] Griewank | Beyond the Oracle: Opportunities of Piecewise Differentiation[END_REF] use settings closer to ours. Piecewise smooth functions, selection derivatives and their variational properties are extensively described in [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]. Our approach differs because we adopt more stringent definitions and rigidity assumptions, which allows in turn for much stronger properties. For instance, we fully treat backward algorithmic differentiation which is the most useful tool in machine learning.

Altogether, our contribution is an accessible and elementary framework for the conservative fields recently introduced in [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF], without explicitly requiring the introduction of semialgebraic geometry and o-minimal structures [START_REF] Van Den Dries | Geometric categories and o-minimal structures[END_REF][START_REF] Coste | An introduction to o-minimal geometry[END_REF].

Stochastic approximation algorithms [START_REF] Robbins | A stochastic approximation method[END_REF][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] are widely used in machine learning contexts [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF][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][START_REF] Lecun | Deep learning[END_REF][START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Chizat | On the global convergence of gradient descent for overparameterized models using optimal transport[END_REF][START_REF] Castera | An Inertial Newton Algorithm for Deep Learning[END_REF]. For example [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF] describes asymptotics of stochastic subgradient algorithms in nonsmooth, nonconvex settings. In contrast, we do not assume access to subgradients and instead explicitly model the behaviour of AD in optimization contexts. Our convergence results are based on [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF], complemented by a new result on "the avoidance of critical traps" in the line of [START_REF] Bianchi | Convergence of constant step stochastic gradient descent for non-smooth non-convex functions[END_REF] in the context of long run convergence.

Notations The ambient space is Euclidean R p . For each k, e k is the k-th vector of the canonical basis. We use D : R m Ñ R q for set valued functions, i.e functions from R m to the subsets of R q . The convex hull of A Ă R p is denoted by convpAq. All proofs are postponed to the Appendix.

Basic piecewise differentiable functions and selection gradient

We introduce a simple but vast class of functions that model rigorously the machine learning models and losses for applications such as deep learning. Examples: The Heaviside function, the index of the largest or k-th largest element in a vector, the sign pattern of a vector in R p which is indexed by integers from 1 to 2 p . Definition 3 (Elementary selection) Let f : R p Þ Ñ R be continuous. We say that f has an elementary (log-exp) selection ps, f 1 , . . . , f m q if s : R p Þ Ñ p1, . . . , mq is an elementary index in I and for i " 1 . . . , m, f i : R p Þ Ñ R are elementary functions in E, such that for all x P R p ,

f pxq " f spxq pxq. (2)
The m `1-uplet ps, f 1 , . . . , f m q is a representation of f , and f admits an elementary (log-exp) selection. The class of such functions is denoted by S log exp or simply here S. This extends to functions from R p to R m by applying a coordinatewise definition with a common elementary index.

Observe that the representation is never unique, both in s and in the sequence f 1 , . . . , f m . The ReLU, hinge loss, maximal entry, k-th largest entry functions are elementary selections. Note also that continuity is part of the definition.

Proposition 1 (Stability of S by ˝, `, ˆ) The class S of elementary selections is stable by composition, sum and product.

The class S is close to the one of piecewise C k functions, see e.g [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF], but it is also much more disciplined since indices and functions are required to satisfy strong "log-exp" rigidity assumptions.

Selection derivative

Functions in S can be associated with a flexible notion of generalized derivative based on the selection structure of the underlying function.

Definition 4 (Selection gradient) (i) Let f : R p Þ Ñ R, in S with selection ps, f 1 , . . . , f m q. We set the selection derivative of f with respect to s to be

p ∇ s f : x Þ Ñ ∇f spxq pxq. (3)
This extends to multivariate outputs by applying the definition coordinatewise, which leads to a notion of a selection Jacobian denoted by p J s . (ii) Given a function f P S, a selection derivative is a derivative of the form (3) for a given representation. In that case a selection derivative of f is merely denoted by p ∇f .

Example: Set for all x P R, f 1 pxq " 0, f 2 pxq " x and spxq " 1 for x ď 0 and spxq " 2 for x ą 0. This this defines the relu function and its selection derivative at 0 is 0. See more in Appendix A.2.

Remark: (a) p ∇f is different from any known notion of subgradient. Set for all x P R, f 1 pxq " 0, f 2 pxq " x and spxq " 1 for x ‰ 0 and sp0q " 2. This defines a elementary selection for the null function however, p ∇ s f p0q " 1. This is the zero function of the introduction. (b) This formalizes what one would obtained by differentiating a code with all decision branches frozen and hence represents the numerical output of AD (see 4). Note that one only needs one branch and do not need to explore all possible outcomes, avoiding combinatorial explosion.

The properties of selection derivatives might seem too liberal at first sight and too disconnected from the original function, but this is not the case as shown below.

Proposition 2 (Integration along segments) Let f : R p Þ Ñ R be in S, with elementary selection ps, f 1 , . . . , f m q. Then f is locally Lipschitz and for all y, x in R p .:

f pyq ´f pxq " ż 1 0 A y ´x, p ∇ s f px `tpy ´xqq E dt Proposition 3 (Gradient almost everywhere) Let f : R p Þ Ñ R be in S, with elementary selection ps, f 1 , . . . , f m q.
There exists sets U 1 , . . . , U N with nonempty interior such that Ť N i"1 clpU i q " R p and for each i " 1, and for all x in the interior of U i , p ∇ s f pxq " ∇f pxq. Furthermore, the U i are solution sets of equations and inequalities involving functions in E.

Remark: Although less transparent, Proposition 2 is not a consequence of Proposition 3. Both results crucially rely on the rigidity of elementary functions in E (Definition 3), not only on their piecewise smoothness. This a central novelty of our approach.

A calculus for selection derivatives

One has an unusual differential calculus: although it does not involve the linearity of some (sub)differential operator, the selection derivative of a sum gives a sum of selection derivatives provided that the selection is refined.

Proposition 4 (Chain rule) Let F : R p1 Þ Ñ R p2 such that each of its coordinate f i , i " 1 . . . p 2 , is in S and g : R p2 Þ Ñ R, g P S. Consider a selection Jacobian for F , p J F : R p1 Þ Ñ R p2ˆp1 x Þ Ñ ¨p ∇f 1 pxq T . . . p ∇f q pxq T ‹ ' (4)
Then g ˝F P S and the function x Þ Ñ p J F pxq T p ∇gpF pxqq is a selection derivative for g ˝F .

Proposition 4 extends readily to the case when the outer function g is multivariate. For example, we have a sum rule p ∇pf `gq " p ∇f `p ∇g for full-domain functions f, g in S. Indeed, if F 1 and F 2 are elementary selections then F 1 ˝F2 P S and p J F1˝F2 " p p J F1 ˝F2 q ˆp J F2 .

(5)

Programs and elementary selections

Numerical programs encode numerical functions by combining elementary functions using a predecessor relation which models program execution. In what follows, m can be seen as an estimate of the memory footprint of a program4 , while p and q the number of inputs and outputs respectively.

Given positive integers m ě p `q, a predecessor relation is a set valued map pr : t1, . . . , mu Ñ t1, . . . , mu such that

• For i P t1, . . . , mu and j P prpiq, j ă i. ' For i P tp `1, . . . , mu, prpiq is nonempty.

A predecessor relation induces a partial order on the set of integers from 1 to m and hence can be represented by a directed acyclic graph [START_REF] Lehman | Mathematics for computer science[END_REF]Theorem 9.4.9]. Given pp, q, mq and a predecessor relation pr, a elementary function sequence G " pg i q m i"p`1 is a set of functions such that g i : R |prpiq| Þ Ñ R, and g i P S, for all i " p `1, . . . , m. A program P is then given by the data P " pp, q, m, pr, Gq , while its evaluation is described in Algorithm 1. We denote by P the set of programs, and P p,q when input-output dimensions have to be made explicit. By definition a program encodes a function, but the representation is not unique. We express this fact below through the canonical surjection F of the introduction.

Algorithm 1: Program evaluation Program data: p, q ě 1, m ě p `q, pr a predecessor relation, G " pg i q m i"p`1 an adapted function sequence.

Input: x " px 1 , . . . x p q 1: for k " p `1, p `2, . . . m do 2:
x k " g k px prpkq q where x prpkq " px i q iPprpkq . 3: end for Return: y :" px j q m j"m´q`1 .

The following proposition illustrates the fact that practitioners implicitly implement selection functions when writing programs.

Proposition 5 (Programs represents elementary selections) Through its input-output correspondence each program P of the form (3) induces a function which is an elementary selection. In other words FpP q P S.

Algorithmic differentiation and a variational model

Algorithmic differentiation is based on the idea of propagating infinitesimal variations in a program P through the chain rule, either forward or backward.

Algorithm 2: Algorithmic differentiation computes selection gradients Program data: p ě 1, m ě p `1, pr a predecessor relation, G " pg i q m i"p`1 an adapted function sequence. Input: variables px 1 , . . . x m q computed by Algorithm 1,

d i " pd i rjsq |prpiq| j"1 " p ∇g i px prpiq q, i " p `1 . . . m 1: Forward mode: 2: Initialize: Bx k Bx " e k , k " 1, . . . , p. 3: for k " p `1, . . . m do 4: Bx k Bx " ÿ jPprpkq Bx j Bx d k rjs
where x " px 1 , . . . , x p q. 5: end for Return: Bxm Bx .

1: Backward mode: 2: Initialize: v " e m 3: for t " m, . . . p `1 do 4:

for j P prptq do 5:

Update coordinate j of v:

vrjs :" vrjs `vrtsd t rjs 6:

end for 7: end for Return: pvr1s, vr2s, . . . , vrpsq.

Consider Algorithm 1, and assume for simplicity that q " 1. The program can be seen as the implementation of m ´p successive transformations on R m , of the form

G k : R m Þ Ñ R m x Þ Ñ x `ek pg k px prpkq q ´xk q,
for k " p `1, . . . , m which belong to S. Algorithm 2 combines gradients dynamically along two modes: forward or backward. Let us describes these two forms.

Fix x P R m . After applying Algorithm 1, for each k, let d k P R m be the selection gradient p ∇g k px prpkq q, appending 0 to non dependant coordinates. A selection Jacobian of G k (at x) is given by

p J G k " I ´ek e T k `ek d T k
Denote by J p P R mˆp , the matrix whose entries are 0, except for diagonal entries which are 1. In Algorithm 2, the forward mode computes

e T m p J Gm . . . p J Gp`1 J p " e T m `I ´em e T m `em d T m ˘. . . `I ´ep`1 e T p`1 `ep`1 d T p`1
˘Jp which is a selection Jacobian thanks to the chain rule in [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF]. On the other hand the backward mode computes J T p `I `dp`1 e T p`1 ˘. . . `I `dm e T m ˘em . This quantity turns out to be the same as the one computed by the forward mode thanks to:

Lemma 1 Let p, m P N, 0 ă p ă m.
where I P R mˆm is the identity matrix and P p P R mˆm denotes the projection on the first p coordinates.

Denote by A : P p,1 Ñ P p,p the algorithmic-differentiation operator. This establishes the following fundamental fact which is at the root of this work. This result asserts that practitioners implicitly implement selection derivatives when writing numerical programs and calling forward or backward AD on these programs.

Theorem 2 (Algorithmic differentiation outputs a selection gradient) Algorithmic differentiation of a given program, i.e., ApP q, outputs a selection derivative of the underlying numerical function. In other words there exists a representation of the numerical function FpP q with elementary index s such that: FpApP qq " p ∇ s FpP q.

5 Algorithmic differentiation at work

Selection derivatives, conservative fields and Clarke subgradient

The asymptotic study of first-order optimization methods implies limiting processes and necessitates thus the introduction of graph closed operators. Given a representation for f , we may construct such a convex-valued mapping pointwise as follows 5 .

Definition 5 (Representation minimal operator) Let f P S with elementary selection ps, f 1 , . . . , f m q. For any x P R p , set Ipxq " ti P t1, . . . , mu , f pxq " f i pxqu. The index closure of p ∇ s f is given by the set valued map

D s f : R p Ñ R p x Ñ conv pt∇f i pxq, i P Ipxquq .
where the double arrows express that the map has values in subsets of R p , much like subgradients, and conv denotes the convex hull.

The role of D s f is to capture all possible outputs of AD including all possible program branches. Of course, due to combinatorial explosion, this quantity is intractable in practice. Its introduction here is only instrumental, we do not use it in algorithms, we just need to access one of its element, for example using a selection derivatives, obtained from AD. A point x satisfying 0 P D s f pxq is called a selection critical point. We will often drop the index s and write D f " D s f . The two following results highlight crucial properties of D f in terms of optimization, they again rely on the rigidity constraint of elementary functions.

Theorem 3 Let f P S with elementary selection ps, f 1 , . . . , f m q and D f be as in Definition 5. Then D f is conservative for f , that is for all absolutely continuous curves γ : r0, 1s Þ Ñ R p , for almost all t P r0, 1s, f ˝γ is differentiable and

d dt f pγptqq " xv, 9 γptqy , @v P D f pγptqq.
The previous result generalizes Proposition 2 by allowing to integrate arbitrary selections along absolutely continuous curves. This connects our work to the general setting of [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF], note that D f has a closed graph thanks to Proposition 6 in Appendix A. [START_REF] Baydin | Automatic differentiation in machine learning: a survey[END_REF].

In [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF], the author considers the essential index set, for each x P R p , S E pxq " ti P t1, . . . , mu , x P clpintpty, f pyq " f i pyquqqu Ă Spxq.

Considering Definition 5 with S E pxq instead of Ipxq leads to the Clarke subgradient, which can also be defined as

B c f pxq " convtd P R p : Dx k P ∆ f , x k Ñ x, ∇f px k q Ñ du
where ∆ f is the dense set of differentiability points of f . While Ipxq can be computed pointwise (check finitely many equalities), it might be very hard to check membership in S E pxq without restrictive qualification conditions on programs [START_REF] Kakade | Provably correct automatic sub-differentiation for qualified programs[END_REF].

Illustration with ReLU and sorting: (a) Set for all x P R, f 1 pxq " 0, f 2 pxq " x, spxq " 1 for x ď 0 and spxq " 2 for x ą 0. This is relu. In this case D f " Brelu, the convex subgradient. (b) Let F : R p Þ Ñ R p to be the sorting function which associates to x a vector P x where P is any permutation such that P x belongs to the set of vectors which values are sorted in descending order coordinatewise. F obviously has an elementary selection and the construction which we have proposed leads to

D F : x Þ Ñ conv tP P ∆, P x " F pxqu ,
where ∆ denotes the set of permutation matrices of size p ˆp. Then D is a conservative mapping for F and it actually corresponds to the Clarke Jacobian.

Convergence of gradient type algorithm and criticality of limit points

Optimization processes in learning are supposed to provide at least a critical point x of the loss, i.e. a point satisfying 0 P B c f pxq. When using AD one enlarges the definition of criticality into 0 P D f pxq and artificial critical points appear, they satisfy 0 R B c f pxq and 0 P D f pxq. Artificial critical points could possibly trap the optimization process in strongly non-optimal situations, we thus have to determine if they have an impact on learning phases.

We consider the problem

min xPR p Jpxq " 1 n n ÿ i"1 f i pxq (7)
where f i : R p Þ Ñ R, f i P S, i " 1, . . . , n. We consider the following algorithm, given x 0 P R p , a sequence of positive step sizes pγ k q kPN and a sequence of iid indices pI k q kPN taken uniformly in the nonempty subsets of t0, . . . , nu,

x k`1 " x k ´γk p ∇f I k px k q where f I " 1 |I| ÿ iPI f i , I Ă t1, . . . , nu. (8)
Note that as discussed in Section 4 selection derivatives can be computed by AD if f i are given by the data of numerical programs as in (3), and could be far from usual notions of subgradients. Hence this algorithm models explicitly the training of a nonsmooth deep network using existing backpropagation implementations. Note that J P S and that 1{n ř n i"1 p ∇f i is a selection gradient for J as stated in Proposition 4, denote by p ∇J this quantity and D J the corresponding set valued field (Definition 5). The following result illustrates that selection critical points are the only attractors for the recursion and that generically such attractors are actually Clarke critical. The first result stands on the theory developed in [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF]. The second parallels developments in [START_REF] Bianchi | Convergence of constant step stochastic gradient descent for non-smooth non-convex functions[END_REF] in the context of long run convergence. The spurious behaviour illustrated in Figure 1 does not affect asymptotics, for typical initialization.

Theorem 4 (Convergence and insignificance of artefacts) Let for all k, γ k " cα k where c P p0, 1s and α k " op1{ log kq and K Ă R p be open. Assume that for all c P p0, 1s and all x 0 P K the sequence in (8) is bounded almost surely.

• For all x 0 P K, almost surely, Jpx k q converges as k tends to infinity and all accumulation points, x, of px k q kPN are selection critical points: 0 P D J pxq.

• For almost all c P p0, 1s, almost all x 0 P K, and almost surely, any accumulation point, x, of px k q kPN is Clarke critical: 0 P B c Jpxq.

Conclusion

The current work departs from existing approaches to nonsmooth algorithmic differentiation in a fundamental way. We propose to study the backward mode of AD, as implemented in machine learning, without any modification. Our theoretical results model thus AD "as is", and our focus is precisely on its unpredictable behavior in a nonsmooth context, addressing an issue which is ubiquitous in machine learning. Our main contribution was to prove that, in a stochastic optimization context, this spurious behavior is essentially harmless from a theoretical point of view, providing justifications for the use of AD outside of its original domain of validity in machine learning.

We achieve our goal by modeling sharply common machine learning functions and their differentiation using selection derivatives, a known concept, which models the way AD differentiates nonsmooth programs. We restrict it to certain classes of elementary functions, opening the possibility to use powerful geometric techniques.

Further questions include convergence rates and complexity issues, hardly tackled at this day, let us mention the attempt of [START_REF] Zhang | On Complexity of Finding Stationary Points of Nonsmooth Nonconvex Functions[END_REF]. Our theory is limited to continuous functions and an interesting venue is to extend it to discontinuous functions, in view of treating ranking operations [START_REF] Blondel | Fast Differentiable Sorting and Ranking[END_REF] ubiquitous in recommendation systems, or more generally differentiating through an argmax [START_REF] Berthet | Learning with differentiable perturbed optimizers[END_REF].

Broader impact

One of the goals of the paper is to raise awareness about an important issue of in the training of ML methods: the spuriousness of AD. To address adequately this issue, we think it is necessary to include algorithmic differentiation explicitly in the study of optimization algorithms, a point of view which is largely ignored by today's machine learning community. This is the appendix for "A mathematical model for automatic differentiation in machine learning".

A A more comprehensive discussion and auxiliary results

A.1 Related work and contribution

The use of backward mode of algorithmic differentiation (AD) for neural network training expanded in the 80's, the most cited reference being [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF]. However the theory applies to much more optimization problems, see for example [START_REF] Griewank | Who invented the reverse mode of differentiation[END_REF]. Indeed, numerical libraries implementing the backward mode of AD were already available in the 90's for FORTRAN code [START_REF] Bischof | ADIFOR-generating derivative codes from Fortran programs[END_REF][START_REF] Bischof | ADIFOR 2.0: Automatic differentiation of Fortran 77 programs[END_REF] or C/C++ code [START_REF] Griewank | Algorithm 755: ADOL-C: a package for the automatic differentiation of algorithms written in C/C++[END_REF], 30 years before the emergence of python libraries. These early implementation could differentiate virtually any code, but their domain of validity, i.e., the setting for which one could predict what the output would be, was restricted to differentiable functions evaluated on their (open) domain of differentiability.

This was well known to the AD community, see for example [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF], and exploring further the domain of validity of AD, beyond mere differentiability, was already a vivid problem.

Let us mention [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] who used notions such as finite selection, "isolated criticalities", stable domain or regular arcs, and argued that "functions given by evaluation procedures are almost everywhere real analytic or stably undefined" where "undefined" meant that a nonsmooth elementary function is used in the evaluation process. For piecewise smooth functions which nonsmoothness can be described using the absolute value function (abs-normal form), [START_REF] Griewank | On stable piecewise linearization and generalized algorithmic differentiation[END_REF] developped a piecewis linearisation formalism and local approximation related to AD, [START_REF] Griewank | On Lipschitz optimization based on gray-box piecewise linearization[END_REF] proposed an AD based bundle type method. These developments are based on the notion of piecewise smooth functions [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF] which we use in this work. More recently, [START_REF] Griewank | Treating artificial neural net training as a nonsmooth global optimization problem[END_REF] applied these techniques to single layer neural network training and [START_REF] Griewank | Beyond the Oracle: Opportunities of Piecewise Differentiation[END_REF] proposed to avoid the usage of subgradient "oracles" in nonsmooth analysis as they are not available in practice. In a similar vein, let us mention [START_REF] Barton | Computationally relevant generalized derivatives: theory, evaluation and applications[END_REF] study lexicographic derivatives, a notion of directional derivatives which satisfy a chain rule making them compatible by forward mode AD, and [START_REF] Zhang | On Complexity of Finding Stationary Points of Nonsmooth Nonconvex Functions[END_REF] who use directional derivatives in the context of local sampling stochastic approximation algorithms for machine learning.

Constraint qualification is known in nonsmooth analysis to ensure favorable behavior of chain rules of differential calculus for nonsmooth objects (see [START_REF] Rockafellar | Variational analysis[END_REF]). These already appeared in the context of piecewise smooth functions of Scholtes with the notion of "essential selections". Such an approach was used in [START_REF] Kakade | Provably correct automatic sub-differentiation for qualified programs[END_REF] to propose an AD algorithm for subgradient computation under constrant qualification. Similarly [START_REF] Griewank | First-and second-order optimality conditions for piecewise smooth objective functions[END_REF] study first and second order optimality, in relation to AD using constraint qualification.

The current work departs from all these approaches in a fundamental way. We propose to study backward mode of AD, as implemented for nonsmooth functions by standard software (e.g. TensorFlow, PyTorch), without any modification, addition of operations or hypotheses. Our theoretical results model AD as implemented in current machine learning libraries. Contrary to previous works, our focus is precisely on the unpredictable behavior of AD in nonsmooth context. Our main contribution is to show that in a stochastic optimization context, this spurious behavior is essentially harmless from a theoretical point of view, providing justifications for the use of AD outside of its original domain of validity in machine learning.

At the time this paper was accepted, we learnt about a paper proposing an analysis close to ours [START_REF] Lee | On Correctness of Automatic Differentiation for Non-Differentiable Functions[END_REF]. The authors show that AD applied to programs involving piecewise analytic continuous functions, under analytic partitions, compute gradients almost everywhere. This is the counterpart of Proposition 3, replacing log-exp elementary function in Definitions 1 and 2, by analytic functions.

A.2 Implementation of relu

The implementation of the relu function used in Figure 1 is given by the function tf.nn.relu in Tensorflow software library [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF]. This implementation corresponds to the selection function described in Section 2 and the same result may be obtained by an explicit implementation of this branching selection as illustrated in the following figure

B Proofs

Proof of Theorem 1: Recall the operator is denoted by B A . Fix a function f , by point (a), the operator

B A f should contain " R p Ñ R p
x Ñ tApP qpxq : FpP q " f, P P Pu Let us show that the graph of the above is R p ˆRp . Assume p " 1 for simplicity. For real numbers r, s, consider the functions f r,s " f `r zerop¨´sq which coincide with f but whose form induces programs P r,s of f . These satisfy FpP r,s q " f and ApP r,s qpsq Q Apf qpsq `r. Since r is arbitrary, B A f psq " R p and since s is arbitrary, we actually have

graph B A f " R p ˆRp .
Since f is arbitrary, we have shown that B A is trivial. l

Proof of Proposition 2:

The proposition is a consequence of Theorem 3 and ( 11) but it admits a more elementary proof which we detail here. Fix x, y P R p . Let us admit the following claim -whose independent proof is given in Section C.

Claim 1 There exists a finite set of numbers 0 " a 0 ă a 1 ă . . . ă a N " 1, such that for all i P 0, . . . N ´1, the function t Þ Ñ spx `tpy ´xqq is constant.

Fix i P 0 . . . , N ´1, and j P 1 . . . m such that f " f j on px `ai py ´xq, x `ai`1 py ´xqq. Since

f j P E p , it

Claim 2

The boundary of each V j has zero measure and cl `Ym i"j intpV j q ˘" R p .

Hence, we may define U 1 , . . . , U N by keeping only those sets with nonempty interior and take their closure. On each set U i , f is identical to f k for some k and the result follows. l

Lemma 2 Let t P I be an elementary index on R p2 and F : R p1 Þ Ñ R p2 with each coordinate in E, then t ˝F is an elementary index on R p1 .

Proof : Fix an arbitrary integer i in the image of t, by Definition 2, there exists elementary functions h 1 , . . . , h J , J P N on R p2 such that tpyq " i if and only if y P K i :" tz P R p2 , h j pzq ˛j 0, j " 1, . . . Ju where ˛j is an equality or inequality sign depending on j. Then tpF pxqq " i if and only if F pxq P K i which is equivalent to say that x P Ki :" tx P R p1 , h j pF pxqq ˛j 0, j " 1, . . . Ju. By Definition 1, h j ˝F is an elementary function for j " 1, . . . , J and i was an arbitrary integer, this shows that we have an elementary index. l

Proof of Proposition 1: Let F : R p1 Þ Ñ R p2 such that each of its coordinate f i , i " 1 . . . p 2 , is in S and g : R p2 Þ Ñ R, g P S. We establish that g ˝F is an elementary selection, the other cases are similar. We may consider all possible intersections of constant index domains across all coordinates of F in t1, . . . , p 2 u. We obtain ps, F 1 , . . . , F m q, an elementary selection for F (each F i : R p1 Þ Ñ R p2 has coordinates in E) . Consider g P S with elementary selection pt, g 1 , . . . , g l q. The composition g ˝F may be written as gpF pxqq " g tpF pxqq pF pxqq " g tpF spxq pxqq pF spxq pxqq.

For each i " 1 . . . , m and j " 1, . . . , l, consider the set

U ij " tx P R p , spxq " i, tpF i pxqq " ju .
Fix pi, jq in t1, . . . , mu ˆt1, . . . , lu, by Lemma 2, t ˝Fi is an elementary index on R p1 . Hence U ij is the solution set of finitely many equalities and inequalities involving functions in E. We associate to the bi-index pi, jq the corresponding set U ij and the function g j pF i pxqq P E. Note that we assumed that the composition is well defined. Identifying each pair pi, jq with a number in t1, . . . , nmu, we obtain an elementary selection for g ˝F and hence g ˝F P S. l

Proof of Proposition 4: The derivation formula follows from the proof argument of Proposition 1, for each pair pi, jq, the function g j ˝Fi is the composition of two C 1 functions and its gradient is given by J Fi ˆ∇g j ˝Fi on U ij . By construction of U ij and definition of the selection derivative, this corresponds to (5) on U ij and the result follows. l

Proof of Lemma 1: We actually prove a slightly stronger result, namely for each i P tp`1, . . . , m1 u

P i `I ´ei`1 e T i`1 `di`1 e T i`1 ˘. . . `I ´em e T m `dm e T m ˘" P i `I `di`1 e T i`1 ˘. . . `I `dm e T m (9)
We argue by exhaustion from i " m ´1 downward to i " p, which is the result of interest. If i " m ´1, we indeed have P m´1 `I ´em e T m `dm e T m ˘" P m´1 `I `dm e T m since P m´1 e m e T m " 0. Now assume that (9) holds true for an index i within tp `1, . . . , m ´1u, then we have

P i´1 `I ´ei e T i `di e T i ˘. . . `I ´em e T m `dm e T m " P i´1 `I ´ei e T i `di e T i ˘`I ´ei`1 e T i`1 `di`1 e T i`1 ˘. . . `I ´em e T m `dm e T m " P i´1 `I ´ei e T i `di e T i ˘Pi `I ´ei`1 e T i`1 `di`1 e T i`1 ˘. . . `I ´em e T m `dm e T m " P i´1 `I `di e T i ˘Pi `I `di`1 e T i`1 ˘. . . `I `dm e T m " P i´1 `I `di e T i ˘`I `di`1 e T i`1 ˘. . . `I `dm e T m ˘,
where step 1 is expanding the product, step 2 is because P i´1 P i " P i´1 and e T i P i " e T i , step 3 combines the fact that P i´1 e i " 0 and (9) which we assumed to be true, the last step uses again the fact that P i´1 P i " P i´1 and e T i P i " e T i . Hence the result holds by exhaustion. l

Proof of Proposition 6: Consider the sequence s k " Spx k q, by taking a subsequence we may assume that s k is constant, say equal to t1, . . . , ru. Hence for all k, v k P conv pt∇f i px k q, i " 1, . . . ruq and f px k q " f i px k q, i " 1, . . . , r. Passing to the limit, we have f pxq " f i pxq, i " 1, . . . , r and hence t1, . . . , ru P Spxq. Furthermore, v P conv pt∇f i pxq, i " 1, . . . ruq Ă D f pxq. l

C o-minimal structures, definability and conservative fields C.1 pR, expq-definability

We recall here the results of geometry that we use in the present work. Some references on this topic are [START_REF] Coste | An introduction to o-minimal geometry[END_REF][START_REF] Van Den Dries | Geometric categories and o-minimal structures[END_REF].

An o-minimal structure on pR, `, ¨q is a collection of sets O " pO p q pPN where each O p is itself a family of subsets of R p , such that for each p P N:

(i) O p is stable by complementation, finite union, finite intersection and contains R p . (ii) if A belongs to O p , then both A ˆR and R ˆA belong to O p`1 ;

(iii) if π : R p`1 Ñ R p is the canonical projection onto R p then, for any A P O p`1 , the set πpAq belongs to O p ; (iv) O p contains the family of real algebraic subsets of R p , that is, every set of the form tx P R p | gpxq " 0u where g : R p Ñ R is a polynomial function; (v) the elements of O 1 are exactly the finite unions of intervals.

A subset of R p which belongs to an o-minimal structure O is said to be definable in O. A function is definable in O whenever its graph is definable in O). A set valued mapping (or a function) is said to be definable in O whenever its graph is definable in O. The terminology tame refers to definability in an o-minimal structure without specifying which structure.

The simplest o-minimal structure is given by the class of real semialgebraic objects. Recall that a set A Ă R p is called semialgebraic if it is a finite union of sets of the form

k č i"1 tx P R p | g i pxq ă 0, h i pxq " 0u
where the functions g i , h i : R p Ñ R are real polynomial functions and k ě 1. The key tool to show that these sets form an o-minimal structure is Tarski-Seidenberg principle which ensures that (iii) holds true.

According to [START_REF] Wilkie | A theorem of the complement and some new o-minimal structures[END_REF], there is an o-minimal structure which contains all semialgebraic sets and the graph of the exponential function, we fix this o-minimal structure and call it O. As a consequence, all functions which can be described by a finite compositional expression involving polynomials, quotients, exponential and logarithms are definable in O. In particular any function f P S is definable in O, which opens the use of powerful geometric tools [START_REF] Coste | An introduction to o-minimal geometry[END_REF][START_REF] Van Den Dries | Geometric categories and o-minimal structures[END_REF] for functions in S. From now on, we call an object definable if it is definable in O.

As detailed in [START_REF] Coste | An introduction to o-minimal geometry[END_REF] the following holds true Proposition 7 (Quantifier elimination) Any first order formula (quantification on variables only) involving definable functions and definable sets describes a definable set.

This allows to prove Claim 1

Proof of Claim 1: The function t Þ Ñ spx `tpy ´xqq is definable and has values in t1, . . . , mu. For each j P t1, . . . , mu, the set S j " tt P r0, 1s, spx `tpy ´xqq " ju is definable, and by (v), it is a finite union of intervals. For each j consider only the endpoints of those intervals with nonempty interior, this provides the desired partition. l

C.2 Properties of definable sets

The tangent space at a point x of a manifold M is denoted by T x M . Given a submanifold 6 M of a finite dimensional Riemannian manifold, it is endowed by the Riemanninan structure inherited from the ambient space. Given f : R p Ñ R and M Ă R p a differentiable submanifold on which f is differentiable, we denote by grad M f its Riemannian gradient or even, when no confusion is possible, grad f . A C r stratification of a (sub)manifold M (of R p ) is a partition S " pM 1 , . . . , M m q of M into C r manifolds having the property that cl M i X M j ‰ H implies that M j is entirely contained in the boundary of M i whenever i ‰ j. Assume that a function f : M Ñ R is given and that M is stratified into manifolds on which f is differentiable. For x in M , we denote by M x the strata containing x and we simply write grad f pxq for the gradient of f with respect to M x .

Stratifications can have many properties, we refer to [START_REF] Van Den Dries | Geometric categories and o-minimal structures[END_REF] and references therein for an account on this question and in particular for more on the idea of a Whitney stratification that we will use repeatedly.

We pertain here to one basic definition: a C r -stratification S " pM i q iPI of a manifold M has the Whitney-(a) property, if for each x P cl M i X M j (with i ‰ j) and for each sequence px k q kPN Ă M i we have:

lim kÑ8 x k " x lim kÑ8 T x k M i " T , / . / - ùñ T x M j Ă T
where the second limit is to be understood in the Grassmanian, i.e., "directional", sense. In the sequel we shall use the term Whitney stratification to refer to a C 1 -stratification with the Whitney-(a) property. The following can be found for example in [START_REF] Van Den Dries | Geometric categories and o-minimal structures[END_REF].

Theorem 5 (Whitney stratification) Let A 1 , . . . , A k be definable subsets of R p , then there exists a definable Whitney stratification pM i q iPI compatible with A 1 , . . . , A k , i.e. such that for each i P I, there is t P t1, . . . ku, such that M i Ă A t .

This allows for example to prove Claim 2

Proof of Claim 2: The sets V 1 , . . . , V m form a definable partition of R p . Consider a Whitney stratification of R p , pM i q iPI compatible with the closure of V 1 , . . . , V m . The boundary of each V i is a finite union of strata of dimension strictly smaller than p and hence has measure zero. The remaining strata (open of maximal dimension) have to be dense in R p since we started with a partition. l C.3 Variational stratification and projection formulas Definition 6 (Variational stratification) Let f : R p Ñ R, be locally Lipschitz continuous, let D : R p Ñ R p be a set valued map and let r ě 1. We say that the couple pf, Dq has a C r variational stratification if there exists a C r Whitney stratification S " pM i q iPI of R p , such that f is C r on each stratum and for all x P R p , Proj T Mx pxq Dpxq " tgrad f pxqu ,

where grad f pxq is the gradient of f restricted to the active strata M x containing x.

The equations [START_REF] Blondel | Fast Differentiable Sorting and Ranking[END_REF] are called projection formulas and are motivated by Corollary 9 in [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF] which states that Clarke subgradients of definable functions have projection formulas.

Let us recall the definition of conservative set-valued mappings from [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF] and one of its characterization. The following combines other results from [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF], where one implication is essentially due to [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF] based on [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF].

Theorem 6 (Characterization of conservativity) Let D : R p Ñ R p be a definable, nonempty compact valued, graph closed set valued field and f : R p Þ Ñ R be a definable locally Lipschitz function.

Then the following are equivalent

• D is conservative for f .

• For any r ě 1, pf, Dq admit a C r variational stratification.

This result allows to prove the following Proof of Theorem 3: We prove that there is a C 1 projection formula (see Theorem 6). For each I Ă t1, . . . , mu, set V I " tx P R p , Spxq " Iu. On each set V I , f pxq " f i pxq for all i P I. These sets are definable, hence, there is a definable Whitney stratification of R p which is compatible with them (Theorem 5). For any C 1 manifold M in the stratification there is an index set I Ă t1, . . . , mu such that for all i P I and all x P M , f pxq " f i pxq and Spxq " I. Since each f i , i P I is C 1 and they agree on M , they represent the same function when restricted to M . Hence they have the same differential on M and since they are all globally C 1 this agrees with the projection of their gradient on the tangent space of M . Hence the projection of D f pxq to the tangent space to M at x is single valued and corresponds to the derivative of f restricted to M . This is sufficient to conclude as this is precisely the variational stratification required by Theorem 6. l

D Convergence to selection critical points

Proof of Theorem 4, first part: We use here the results on conservative fields developed in [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF].

To prove the theorem it suffices to establish that:

• D J is a conservative field for J • the number of D J critical values are finite.

The first point is Theorem 6 while the second one is the consequence of the latter and the definability of the couple f, D f , see Proposition 8 (ii). To conclude it suffices to apply the convergence results in [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF]Theorem 9]. l

Proof of Theorem 4, second part: This result is a consequence of the more general Theorem 7 established in Section E. Let F be the finite set given in Theorem 7, the set tc P p0, 1s, Dk P N, cγ k P F u, is countable, and hence has zero measure. So for almost all c P p0, 1s, tcγ k u kPN does not intersect F . Using Theorem 7, there is a zero measure set N such that any initialization outside N provides almost surely a subgradient sequence. By hypothesis, for almost every x 0 P KzN , the sequence is bounded almost surely and the result follows from Theorem 7. l

E Artificial critical points

Being given a Lipschitz continuous function on R p and a conservative field D, one has two types of D-critical points:

• Clarke critical points: B c f pxq Q 0, we denote the set of these points by crit c f • Artificial critical points B c f pxq S 0 and Dpxq Q 0, we denote this set by crit a f Critical values are defined accordingly as images of critical points. Proof : By definability of crit a f , the number of connected components is finite.

If C i had full dimension it would contain a non trivial ball on which f should be constant by the integral property. This would in turn imply that the points in the ball would also be local minimum and thus Clarke critical, which is impossible.

To see that the critical values are in finite number it suffices to evoke the fact that Clarke critical values are finite [START_REF] Bolte | Clarke subgradients of stratifiable functions[END_REF] and use that artificial critical values are in finite number.

By definability the connected components are arcwise-connected with piecewise C 1 paths. Using the integral property this shows f is constant on C i .

(iii) is obvious since local minimum or maximum are Clarke critical. l

As explained in the introduction, artificial critical points are "computing artefacts", whence their names. For algorithmic differentiation the "gradient" provided by a program is zero while the point might even be a smooth non critical point. We consider the setting of the mini-batch algorithm of the last section.

Theorem 7 Assume that each f 1 , . . . , f n belongs to S. There exists a finite subset of steps F Ă p0, `8q and a zero measure meager subset N of R p , such that for any positive sequence γ k " op1{ log kq avoiding values in F , and any almost surely bounded sequence with initial condition in R p zN , we have

• Jpx k q converges towards a Clarke critical value almost surely,

• the cluster points of x k are Clarke critical point almost surely.

Proof : The proof is twofold. We first prove that the set of initial conditions leading to an artificial critical point or more generally to a non differentiability point within a finite time is "small". We then use this fact to conclude.

Claim 3 Let g : R p Ñ R be a definable differentiable function. Set, for λ ą 0, Φ λ " λId ´∇g,

where Id denotes the identity. There exists a finite set F in p0, `8q such that, @λ P p0, `8qzF, @Z Ă R p definable , dim Z ă p ñ dim Φ ´1 λ pZq ă p.

Proof of the claim. Denote by L the set of points where g is twice differentiable so that L is dense and definable. Denote by λ 1 , . . . , λ p : L Ñ R a representation of the eigenvalues of ∇ 2 g. Refine L to be contained in the common domain of differentiability for each λ i , L remains open and dense. By the definable Sard's theorem the critical values of each function λ i is finite, so that the set of all these values which we denote by F is itself finite.

Take a positive real λ R F and consider the set K λ :" tx P L : Φ 1 λ pxq " λId ´∇2 gpxq is not invertibleu.

By diagonalization, we see that the determinant of Φ 1 λ pxq is

d ź i"1
pλ ´λi pxqq for any x, thence

K λ Ă m ď i"1 tx P L, λ i pxq " λu.
Since λ is a regular value for each λ i the previous set is a finite union of manifolds of dimension p ´1, see e.g., [START_REF] Coste | An introduction to o-minimal geometry[END_REF]. This implies that the set R p zK λ " tx P L : Φ 1 λ pxq is invertible u is dense. Using the above, we deduce that there exists finitely many open connected subsets U 1 , . . . , U r Ă L of R p zK λ such that U 1 Y . . . Y U r is dense in L and thus in R p . Take now Z Ă R p definable with dim Z ă p. Assume towards a contradiction that there exists a nonempty open ball B in Φ ´1 λ pZq. In that case B must have a nonempty intersection with some U i0 . The set Φ λ pB X U i0 q is open because Φ λ is a diffeomorphism on U i on its image. Since on the other hand we have Φ λ pB X U i0 q Ă Z, we have a contradiction and the claim is proved. l Set Φ k I,j " Id ´γk ∇f I,j , with I Ă t1, . . . , mu, j P t1, . . . , m I u and Id denotes the identity. Applying Claim 3, we can find a finite set F for which γ k R F implies that each Φ k I,j has the property [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF]. Indeed, for each I Ă t1, . . . , mu, j P t1, . . . , m I u, there is F I,j Ă R finite such that f I,j has property [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF]. Since the subsets I are in finite number and each m I is finite, the set F " Ť IĂt1,...,mu Ť jPt1,...,m I u F I,j , is also finite. For each k P N, I Ă t1, . . . , mu, j P t1, . . . , m I u. Remark that if γ k R F then Φ k I,j has property [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF]. For k ď k 0 fixed, let us consider the finite set of definable mappings defined by

Ψ k0 :" # k ź j"1 Φ j
Ij ,ij : k ď k 0 , I j Ă t1, . . . , nu, i j P t1, . . . , m Ij u + .

We now assume that γ k R F, @k ě 0, so that each mapping in Ψ k0 has the property [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF] and N k0 :" tx P R p : Dk ď k 0 , DΦ P Ψ k pxq P C Y Singu

These are initial conditions in U leading to an artificial critical or a non-differentiability point within U before time k 0 .

We can also write

N k0 Ă ď ΦPΨ k 0 Φ ´1 pC Y Singq .
From stratification arguments we know that Sing has a dimension lower than p ´1. On the other hand, C has dimension strictly lower than p by Proposition 3. Claim 3 applies and yields dim Φ ´1 pC Y Singq ă p for all Φ P Φ k0 . As a consequence N k0 is closed with nonempty interior and so does N :" Y kPN N k by Baire's theorem. Similarly N has zero measure as a countable union of zero measure sets.

This proves that any sequence with initial condition out of N must remain in the zone of differentiability of J as well as all f I . In particular if I is taken uniformly at random among possible subsets, for all x R N , we have E I r p ∇f I pxqs " p ∇Jpxq " ∇Jpxq " B c Jpxq, so that these specific sequences can also be seen as stochastic subgradient sequences for J. To be more specific, the sequence x k can be seen as one of the sequence generated by the algorithm

y k`1 P y k ´γk B c Jpy k q ` k
where k is a random noise with zero mean. Using general results [START_REF] Davis | Stochastic subgradient method converges on tame functions[END_REF][START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF], we know that y k sequences, when bounded almost surely, have limit points which are Clarke critical. l

34th

  Conference on Neural Information Processing Systems (NeurIPS 2020), Vancouver, Canada.

Figure 1 :

 1 Figure 1: Top: AD applied to relu and two different implementations of the same function. Bottom: Algorithmic differentiation of a constant function, creation of artificial critical point or arbitrary derivatives at prescribed arguments for the sine function.

Figure 2 :

 2 Figure 2: Left: Algorithmic differentiation applied to programs combining smooth functions in a smooth way, the diagram commutes. Right: Algorithmic differentiation in nonsmooth settings, connection with known notion of generalized derivative is much less clear.

  One can imagine an equivalent implementation of relu with a slightly different branching involving a strict inequality, that would correspond to an equivalent implementation of the same function, but the computed derivative at 0 is different due to the implementationA.3 Auxiliary results and remarksRemark 1 (Elementary piecewise differentiable functions) (a) The building blocks in the construction of S in Definition 3 could be modified and adapted to other needs. Besides, the results we present in this article would remain true if we added real analytic functions restricted to compact sets. (b) Note also that in Definition 3, functions are actually real analytic on their (open) domain of definition. Yet their extension might not be analytic, as for instance the function f : x ‰ 0 Ñ expp´1{x 2 q extended by f p0q " 0. (c) The construction of elementary piecewise functions in Definition 3, does not coincide in general with some natural minimal o-minimal, but are contained in a larger such structure. For instance, when the basic bricks are polynomial functions, we obtain the field of rational functions which differs from the set of semi-algebraic functions.Proposition 6 (D f has a closed graph) As k Ñ 8, assume that x k Ñ x P R p and v k P D f px k q, v k Ñ v. Then v P Dpxq.

For

  each I Ă t1, . . . , nu, we denote by f I,1 , . . . , f I,m I the bricks attached to f I where m I ě 1. Denote by Sing the set of points on which at least one f I is non differentiable and C the set of points for which p ∇f I ‰ ∇f I for at least one I. By Proposition 3 and definability, Sing and C are finite unions of manifolds of dimension at most p ´1.

  Examples include polynomials, logistic loss, boosting loss, Gaussian likelihood. Observe that the corresponding functions are C 8 smooth on their open domains. Note also that if log and exp are not present we obtain the field of rational functions. See Remark 1 in Appendix A.3.Definition 2 (Elementary index) s : R p Þ Ñ t1, . . . , mu is an elementary (log-exp) index if the set tx P R p , spxq " iu is the solution set of a finite number of inequalities and equalities involving elementary functions on R p . The set of such functions is denoted by I (for any input dimensions p).

	Definition 1 (Elementary (log-exp) functions) Elementary (log-exp) functions are functions on R p
	described by a finite compositional expression involving basic operations, `, ´, ˆ, { as well as affine
	mappings, exponential and logarithms, inside their domain of definition. We denote by E the set of
	elementary functions in any dimension p.

  Assume that for i " p `1, . . . , m we have d i P R m . Then we have P

p `I ´ep`1 e T p`1 `dp`1 e T p`1 ˘. . . `I ´em e T m `dm e T m ˘" P p `I `dp`1 e T p`1 ˘. . . `I `dm e T m

  is C 1 and we have by the fundamental theorem of integral calculus f px `ai`1 py ´xqq ´f px `ai py ´xqq "

				ż ai`1
				x∇f j px `tpy ´xqq, y ´xy dt
				ai
				ż ai`1	A
				"	p ∇f px `tpy ´xqq, y	´xE dt.
				ai
	The conclusion follows because
	N ´1
	ÿ
	f pyq ´f pxq "		f px `ai`1 py ´xqq ´f px `ai py ´xqq
	i"0
	N ´1 ÿ	ż ai`1	A
	"			p ∇f px `tpy ´xqq, y	´xE dt
	i"0	ai
	ż 1	A
	"		p ∇f px `tpy ´xqq, y	´xE dt.
	0

l Proof of Proposition 3: Constructs the sets U i by considering sets V j " tx P R p , spxq " ju, j 1 . . . m, the proof of the following claim is postponed to Section C.

  Definition 7 (Conservative set-valued mappings) Let f be a Lipschitz continuous function. A set valued vector field D is called conservative if for any absolutely continuous path γ : r0, 1s Þ Ñ R, we have Equivalently D is conservative for f , if for all absolutely continuous curves γ : r0, 1s Þ Ñ R p , for almost all t P r0, 1s, f ˝γ is differentiable and

			ż 1				ż 1
	f pγp1qq ´f pγp0qq "	0	min vPDpγptqq	xv, 9 γptqy dt "	0	max vPDpγptqq	xv, 9 γptqy dt.	(11)
	d dt	f pγptqq " xv, 9 γptqy ,	@v P Dpγptqq.

In the usual mathematical sense.

We consider programs which do not overwrite values in memory

Minimality relates to the representation of the function, not the function itself. This is the minimal convex-valued operator, constructed pointwise and guaranteed to be graph-closed.

We only consider embedded submanifolds

Acknowledgments and Disclosure of Funding

The authors acknowledge the support of ANR-3IA Artificial and Natural Intelligence Toulouse Institute, Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbers FA9550-19-1-7026, FA9550-18-1-0226, and ANR MasDol. J. Bolte acknowledges the support of ANR Chess, grant ANR-17-EURE-0010 and ANR OMS. The authors would like to thank anonymous referees for careful reading of this work and useful suggestions. The authors would like to thank N. Asher and S. Gerchinovitz for useful discussions. We also warmly thank J. Malick who triggered this research.