A mathematical model for automatic differentiation in machine learning

Introduction

Optimization algorithms based on backpropagation oracles, and more generally on automatic or algorithmic differentiation (AD) , are one of the most widely used training tools for modern learning architectures . They often rely on popular numerical implementations as TensorFlow or PyTorch . However, for nonsmooth, nonconvex losses, AD does not have a stable theory , matching the actual practice. We wish to present a simple mathematical framework addressing this issue. Let us progressively explain our approach.

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 functionIn the usual mathematical sense. of its arguments.

It is tempting to identify both, but functions can be represented by different programs. This defines a surjection $\mathcal{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, $\mathcal{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 $\mathcal{F}$. Remarkably, if we restrict ourselves to programs $\mathcal{P}$ which only smoothly combine smooth functions, then we have the following fundamental relation, depicted in Figure 2:

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 $\partial^{A}$, should:

encompass the outputs of algorithmic differentation for all functions

Unfortunately such an operator does not exist:

There is no nontrivial operator on functions satisfying $(a)$ and $(b)$.

2 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.

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 . 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 . 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 , 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 . In particular, considers qualification conditions allowing to construct AD algorithms which compute proper Clarke subgradients . However qualification is extremely hard to check and almost impossible to enforce in practice. Let us also mention 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.

use settings closer to ours. Piecewise smooth functions, selection derivatives and their variational properties are extensively described in . 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 , without explicitly requiring the introduction of semialgebraic geometry and o-minimal structures .

Stochastic approximation algorithms are widely used in machine learning contexts . For example 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 , complemented by a new result on “the avoidance of critical traps” in the line of in the context of long run convergence.

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 include polynomials, logistic loss, boosting loss, Gaussian likelihood. Observe that the corresponding functions are $C^{\infty}$ 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.

Observe that the representation is never unique, both in $s$ and in the sequence $f_{1},\ldots,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.

\circ,+,\times) The class $\mathcal{S}$ of elementary selections is stable by composition, sum and product.

The class $\mathcal{S}$ is close to the one of piecewise $C^{k}$ functions, see e.g , but it is also much more disciplined since indices and functions are required to satisfy strong “log-exp” rigidity assumptions.

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

This extends to multivariate outputs by applying the definition coordinatewise, which leads to a notion of a selection Jacobian denoted by $\widehat{J}^{s}$. (ii) Given a function $f\in\mathcal{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 $\widehat{\nabla}f$.

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.

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

2 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.

Then $g\circ F\in\mathcal{S}$ and the function $x\mapsto\widehat{J}_{F}(x)^{T}\widehat{\nabla}g(F(x))$ is a selection derivative for $g\circ F$.

Proposition 4 extends readily to the case when the outer function $g$ is multivariate. For example, we have a sum rule $\widehat{\nabla}(f+g)=\widehat{\nabla}f+\widehat{\nabla}g$ for full-domain functions $f,g$ in $\mathcal{S}$. Indeed, if $F_{1}$ and $F_{2}$ are elementary selections then $F_{1}\circ F_{2}\in\mathcal{S}$ and

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 programWe consider programs which do not overwrite values in memory, while $p$ and $q$ the number of inputs and outputs respectively.

By definition a program encodes a function, but the representation is not unique. We express this fact below through the canonical surjection $\mathcal{F}$ of the introduction.

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

Through its input-output correspondence each program $P$ of the form (3) induces a function which is an elementary selection. In other words $\mathcal{F}(P)\in\mathcal{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.

for $k=p+1,\ldots,m$ which belong to $\mathcal{S}$. Algorithm 2 combines gradients dynamically along two modes: forward or backward. Let us describes these two forms.

which is a selection Jacobian thanks to the chain rule in (5). On the other hand the backward mode computes

This quantity turns out to be the same as the one computed by the forward mode thanks to:

Denote by $\mathcal{A}:\mathcal{P}_{p,1}\to\mathcal{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.

Algorithmic differentiation of a given program, i.e., $\mathcal{A}(P)$, outputs a selection derivative of the underlying numerical function. In other words there exists a representation of the numerical function $\mathcal{F}(P)$ with elementary index $s$ such that:

Algorithmic differentiation at work

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 followsMinimality 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..

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\in D^{s}_{f}(x)$ 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.

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 , note that $D_{f}$ has a closed graph thanks to Proposition 6 in Appendix A.3.

Considering Definition 5 with $S_{E}(x)$ instead of $I(x)$ leads to the Clarke subgradient, which can also be defined as

where $\Delta_{f}$ is the dense set of differentiability points of $f$. While $I(x)$ can be computed pointwise (check finitely many equalities), it might be very hard to check membership in $S_{E}(x)$ without restrictive qualification conditions on programs .

where $\Delta$ denotes the set of permutation matrices of size $p\times p$. Then $D$ is a conservative mapping for $F$ and it actually corresponds to the Clarke Jacobian.

2 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\in\partial^{c}f(x)$. When using AD one enlarges the definition of criticality into $0\in D_{f}(x)$ and artificial critical points appear, they satisfy $0\not\in\partial^{c}f(x)$ and $0\in D_{f}(x)$. 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.

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\in\mathcal{S}$ and that $1/n\sum_{i=1}^{n}\widehat{\nabla}f_{i}$ is a selection gradient for $J$ as stated in Proposition 4, denote by $\widehat{\nabla}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 . The second parallels developments in in the context of long run convergence. The spurious behaviour illustrated in Figure 1 does not affect asymptotics, for typical initialization.

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 . Our theory is limited to continuous functions and an interesting venue is to extend it to discontinuous functions, in view of treating ranking operations ubiquitous in recommendation systems, or more generally differentiating through an argmax .

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.

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.

References

Appendix A A more comprehensive discussion and auxiliary results

The use of backward mode of algorithmic differentiation (AD) for neural network training expanded in the 80’s, the most cited reference being . However the theory applies to much more optimization problems, see for example . Indeed, numerical libraries implementing the backward mode of AD were already available in the 90’s for FORTRAN code or C/C++ code , 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 , and exploring further the domain of validity of AD, beyond mere differentiability, was already a vivid problem.

Let us mention 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), developped a piecewis linearisation formalism and local approximation related to AD, proposed an AD based bundle type method. These developments are based on the notion of piecewise smooth functions which we use in this work. More recently, applied these techniques to single layer neural network training and 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 study lexicographic derivatives, a notion of directional derivatives which satisfy a chain rule making them compatible by forward mode AD, and 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 ). These already appeared in the context of piecewise smooth functions of Scholtes with the notion of “essential selections”. Such an approach was used in to propose an AD algorithm for subgradient computation under constrant qualification. Similarly 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 . 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.

The implementation of the $relu$ function used in Figure 1 is given by the function tf.nn.relu in Tensorflow software library . 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

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 is different due to the implementation

A.3 Auxiliary results and remarks

(a) The building blocks in the construction of $\mathcal{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\neq 0\to\exp(-1/x^{2})$ extended by $f(0)=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.

Appendix B Proofs

Proof of Theorem 1: Recall the operator is denoted by $\partial^{A}$. Fix a function $f$, by point (a), the operator $\partial^{A}f$ should contain

Since $f$ is arbitrary, we have shown that $\partial^{A}$ is trivial. $\Box$

There exists a finite set of numbers $0=a_{0}<a_{1}<\ldots<a_{N}=1$, such that for all $i\in 0,\ldots N-1$, the function $t\mapsto s(x+t(y-x))$ is constant.

Fix $i\in 0\ldots,N-1$, and $j\in 1\ldots m$ such that $f=f_{j}$ on $(x+a_{i}(y-x),x+a_{i+1}(y-x))$. Since $f_{j}\in\mathcal{E}_{p}$, it is $C^{1}$ and we have by the fundamental theorem of integral calculus

Hence, we may define $U_{1},\ldots,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. $\Box$

For each $i=1\ldots,m$ and $j=1,\ldots,l$, consider the set

Proof of Proposition 4: The derivation formula follows from the proof argument of Proposition 1, for each pair $(i,j)$, the function $g_{j}\circ F_{i}$ is the composition of two $C^{1}$ functions and its gradient is given by $J_{F_{i}}\times\nabla g_{j}\circ F_{i}$ 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. $\Box$

Proof of Lemma 1: We actually prove a slightly stronger result, namely for each $i\in\{p+1,\ldots,m-1\}$

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

since $P_{m-1}e_{m}e_{m}^{T}=0$. Now assume that (9) holds true for an index $i$ within $\{p+1,\ldots,m-1\}$, then we have

where step 1 is expanding the product, step 2 is because $P_{i-1}P_{i}=P_{i-1}$ and $e_{i}^{T}P_{i}=e_{i}^{T}$, 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_{i}^{T}P_{i}=e_{i}^{T}$. Hence the result holds by exhaustion. $\Box$

Appendix C o-minimal structures, definability and conservative fields

We recall here the results of geometry that we use in the present work. Some references on this topic are .

the elements of $\mathcal{O}_{1}$ are exactly the finite unions of intervals.

According to , 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 $\mathcal{O}$. As a consequence, all functions which can be described by a finite compositional expression involving polynomials, quotients, exponential and logarithms are definable in $\mathcal{O}$. In particular any function $f\in\mathcal{S}$ is definable in $\mathcal{O}$, which opens the use of powerful geometric tools for functions in $\mathcal{S}$. From now on, we call an object definable if it is definable in $\mathcal{O}$.

Any first order formula (quantification on variables only) involving definable functions and definable sets describes a definable set.

Proof of Claim 1: The function $t\mapsto s(x+t(y-x))$ is definable and has values in $\{1,\ldots,m\}$. For each $j\in\{1,\ldots,m\}$, the set $S_{j}=\{t\in,\,s(x+t(y-x))=j\}$ 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. $\Box$

C.2 Properties of definable sets

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 .

C.3 Variational stratification and projection formulas

The equations (10) are called projection formulas and are motivated by Corollary 9 in which states that Clarke subgradients of definable functions have projection formulas.

Let us recall the definition of conservative set-valued mappings from and one of its characterization.

The following combines other results from , where one implication is essentially due to based on .

For any $r\geq 1$, $(f,D)$ admit a $C^{r}$ variational stratification.

This result allows to prove the following

Appendix D Convergence to selection critical points

Proof of Theorem 4, first part: We use here the results on conservative fields developed in . To prove the theorem it suffices to establish that:

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 [12, Theorem 9]. $\Box$

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

Appendix E Artificial critical points

Critical values are defined accordingly as images of critical points.

$f(C_{i})$ is a singleton, and as a consequence the $D$ critical values of $f$ are in finite number,

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 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. $\Box$

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.

$J(x^{k})$ 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.

where $Id$ denotes the identity. There exists a finite set $F$ in $(0,+\infty)$ such that,

Take a positive real $\lambda\notin F$ and consider the set

By diagonalization, we see that the determinant of $\Phi^{\prime}_{\lambda}(x)$ is $\displaystyle\prod_{i=1}^{d}(\lambda-\lambda_{i}(x))$ for any $x$, thence

For $k\leq k_{0}$ fixed, let us consider the finite set of definable mappings defined by

We now assume that $\gamma_{k}\notin F,\forall k\geq 0,$ so that each mapping in $\Psi_{k_{0}}$ has the property $\eqref{e:dim}$ and

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

where $\epsilon_{k}$ is a random noise with zero mean. Using general results , we know that $y_{k}$ sequences, when bounded almost surely, have limit points which are Clarke critical. $\Box$