Aloïs Brunel

Damiano Mazza

Michele Pagani

France Deepomatic

Damiano Mazza

Backpropagation in the Simply Typed Lambda-Calculus with Linear Negation

Keywords: CCS Concepts:, Theory of computation → Program semantics, Theory and algorithms for application domains, Differentiable Programming, Lambda-calculus, Linear Logic

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

INTRODUCTION

In the past decade there has been a surge of interest in so-called deep learning, a class of machine learning methods based on multi-layer neural networks. The term "deep" has no formal meaning, it is essentially a synonym of "multi-layer", which refers to the fact that the networks have, together with their input and output layers, at least one internal (or "hidden") layer of artificial neurons. These are the basic building blocks of neural networks: technically, they are just functions R m → R of the form (x 1 , . . . , x m ) → σ m i=1 w i • x i , where σ : R → R is some activation function and w 1 , . . . , w m ∈ R are the weights of the neuron. A layer consists of several unconnected artificial neurons, in parallel.

The simplest way of arranging layers is to connect them in cascade, every output of each layer feeding into every neuron of the next layer, forming a so-called feed-forward architecture. Feed-forward, multi-layer neural networks are known to be universal approximators: any continuous function f : K → R with K ⊆ R k compact may be approximated to an arbitrary degree of precision by a feed-forward neural network with one hidden layer, as long as the weights are set properly [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF]. This leads us to the question of how to efficiently train a neural network, i.e., how to find the right weights as quickly as possible.

Albeit in theory the simplest architecture suffices for all purposes, in practice more complex architectures may perform better, and different architectures may be better suited for different purposes. We thus get to another basic question of deep learning: how to select and, if necessary, modify or design network architectures adapted to a given task.

Since the quality of an architecture is also judged in terms of training efficiency, this problem is actually interlaced with the previous one.

The first question is generally answered in terms of the gradient descent algorithm (or some variant thereof). This finds local minima of a function G : R n → R using its gradient ∇G, i.e., the vector of the partial derivatives of G (Equation 1). The algorithm starts by choosing a point w 0 ∈ R n . Under certain assumptions, if ∇G(w 0 ) is close to zero Authors' addresses: Aloïs Brunel, Deepomatic, France, alois.brunel@gmail.com; Damiano Mazza, CNRS, UMR 7030, LIPN, Université Sorbonne Paris Nord, France, Damiano.Mazza@lipn.univ-paris13.fr; Michele Pagani, IRIF UMR 8243, Université de Paris, CNRS, France, pagani@irif.fr.

Manuscript submitted to ACM 1

x 1 x 2 z 1 z 2 y - • sin let z 1 = x 1 -x 2 in let z 2 = z 1 • z 1 in sin z 2
Fig. 1. A computational graph with inputs x 1 , x 2 and output y, and its corresponding term. Nodes are drawn as circles, hyperedges as triangles. The output y does not appear in the term: it corresponds to its root.

then w 0 is within a sufficiently small neighborhood of a local minimum. Otherwise, we know that G decreases most sharply in the opposite direction of ∇G(w 0 ), and so the algorithm sets w 1 := w 0 -ρ∇G(w 0 ) for a suitable step rate ρ > 0, and repeats the procedure from w 1 . In the case of deep learning, a neural network (with one output) induces a function h(w, x) : R n+k → R, where k is the number of inputs and n the number of weights of all the neurons in the network. By fixing w ∈ R n and making x vary over a set of training inputs, we may measure how much h differs from the target function f : R k → R when the weights are set to w. This defines an error function G : R n → R, the minimization of which gives a way of finding the desired value for the weights. The training process thus becomes the iteration, over and over, of a gradient computation, starting from some initial set of weights w 0 .

So, regardless of the architecture, efficiently training a neural network involves efficiently computing gradients.

The interest of gradient descent, however, goes well beyond deep learning, into fields such as physical modeling and engineering design optimization, each with numerous applications. It is thus no wonder that a whole research field, known as automatic differentiation (AD for short), developed around the computation of gradients and, more generally, Jacobians 1 [START_REF] Güneş Baydin | Automatic Differentiation in Machine Learning: a Survey[END_REF]. In AD, the setting of neural networks is generalized to computational graphs, which are arbitrarily complex compositions of nodes computing basic functions and in which the output of a node may be shared as the input of an unbounded number of nodes. Fig. 1 gives a pictorial representation of a simple computational graph made of only three basic functions (subtraction, multiplication and sine), with inputs x 1 and x 2 . Notice that the computation of x 1 -x 2 is shared by the two factors of the multiplication. Neural networks are special cases of computational graphs.

The key idea of AD is to compute the gradient of a computational graph by accumulating in a suitable way the partial derivatives of the basic functions composing the graph. This rests on the mathematical principle known as chain rule, giving the derivative of a composition f • д from the derivatives of its components f and д (Equation 3). There are two main "modes" of applying this rule in AD, either forward, propagating derivatives from inputs to outputs, or backwards, propagating derivatives from outputs to inputs. As will be explained in Sect. 2, if G is a computational graph becomes astronomical (millions, or even billions) while m = 1, hence the reverse mode is the method of choice and specializes in what is called the backpropagation algorithm (Sect. 2.4), which has been a staple of deep learning for decades [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF]]. Today, AD techniques are routinely used in the industry through deep learning frameworks such as TensorFlow [START_REF] Abadi | TensorFlow: A System for Large-Scale Machine Learning[END_REF] and PyTorch [START_REF] Paszke | Automatic differentiation in PyTorch[END_REF].

The interest of the programming languages (PL) community in AD stems from the second deep learning question mentioned above, namely the design and manipulation of (complex) neural network architectures. As it turns out, 1 The generalization of the gradient to the case of functions R n → R m with m > 1.

these are being expressed more and more commonly in terms of actual programs, with branching, recursive calls or even higher-order primitives, like list combinators such as map or fold, to the point of yielding what some call a generalization of deep learning, branded differentiable programming [START_REF] Lecun | Deep Learning est mort[END_REF]]. Although these programs always reduce, in the end, to computational graphs, these latter are inherently static and therefore inadequate to properly describe something which is, in reality, a dynamically-generated neural network. Similarly, traditional AD falls short of providing a fully general foundation for differentiable programming because, in order to compute the gradient of an a priori arbitrary (higher order) program, it forces us to reduce it to a computational graph first. In PL-theoretic terms, this amounts to fixing a reduction strategy, which cannot always be optimal in terms of efficiency. There is also a question of modularity: if gradients may only be computed by running programs, then we are implicitly rejecting the possibility of computing gradients modularly, because a minor change in the code might result in having to recompute everything from scratch, which is clearly unsatisfactory.

This paper is a contribution to the theoretical foundations of differentiable programming. We define a compositional program transformation ← -D (Table 3) extending the backpropagation algorithm from computational graphs to general simply typed λ-terms. Our framework is purely logical and therefore offers the possibility of importing tools from semantics, type systems and rewriting theory. The benefit is at least twofold, in the form of (1) a soundness proof (Theorem 5.6 and Corollary 5.7), which relies on the logical/compositional definition of

← -D

and the semantics;

(2) a complexity analysis, which hinges on rewriting in the form of the linear substitution calculus [START_REF] Accattoli | An Abstract Factorization Theorem for Explicit Substitutions[END_REF] and which guarantees that generalized backpropagation is at least as efficient as the standard algorithm on computational graphs.

Although the soundness proof is based on a fixed strategy (reducing to a computational graph first and then computing the gradient), the confluence of our calculus guarantees us that the gradient may be computed using any other reduction strategy, thus allowing arbitrary evaluation mechanisms for executing backpropagation in a purely functional setting, without necessarily passing through computational graphs. Sect. 6 discusses the benefits of this in terms of efficiency.

On a similar note, compositionality ensures modularity: for instance, if p = tu, i.e., program p is composed of subprograms t and u, then

← - D (p) = ← - D (t) ← - D (
u) and ← -D (t) and ← -D (u) may be computed independently, so that if p is modified into tu ′ , the computation of ← -D (t) may be reused. Modularity also opens the way to parallelization, another potential avenue to efficiency.

Finally, let us stress that the transformation ← -D is remarkably simple: on "pure" λ-terms, i.e., not containing any function symbol corresponding to the basic nodes of computational graphs, ← -D is the identity, modulo a change in the types. In particular, ← -D maps a datatype constructor/destructor (cons, head, tail, etc.) or a typical higher order combinator (map, fold, etc.) to itself, which makes the transformation particularly easy to compute and analyze. We see this as a theoretical counterpart to the use of operator overloading in implementations of AD.

From a more abstract perspective, all these properties (including compositionality) may be succinctly summarized in the fact that ← -D is a cartesian closed 2-functor or, better, a morphism of cartesian 2-multicategories, obtained by freely lifting to λ-terms a morphism defined on computational graphs. However, such a viewpoint will not be developed in this paper.

Related work. Our main source of inspiration is [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF], where a program transformation

← - D is proposed
as a compositional extension of symbolic backpropagation to functional programs. We summarize their approach in Sect. 3, let us concentrate on the main differences here:

3

(1) the transformation ← -D uses references and delimited continuations, while our transformation ← -D is purely functional and only relies on a linear negation primitive on the ground type. Albeit [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]] do mention that a purely functional version of their transformation may be obtained by encoding the memory inside the type of the continuation, this encoding adds a sequentialization (introduced by the order of the memory updates)

which is absent in ← -D and which makes our transformation more amenable to parallelization.

(2) The transformation ← -D applies to a Turing-complete programming language, while we focus on the much more restrictive simply-typed λ-calculus. However, no soundness proof for ← -D is given, only testing on examples. This brings to light a difference in the general spirit of the two approaches: the paper [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]] is mainly experimental, it comes with a deep learning framework (Lantern) and with a case study showing the relevance of this line of research. Our approach is complementary because it is mainly theoretical: we focus on proving the soundness and complexity bound of a non-trivial differentiable programming framework, "non-trivial" in the sense that it is the simplest exhibiting the difficulty of the task at hand. To the best of our knowledge, such foundational questions have been completely neglected, even for "toy" languages, and our paper is the first to address them.

The earliest work performing AD in a functional setting is [START_REF] Barak | Reverse-mode AD in a Functional Framework: Lambda the Ultimate Backpropagator[END_REF]. Their motivations are broader than ours: they want to define a programming language with the ability to perform AD on its own programs.

To this end, they endow Scheme with a combinator ← -J computing the Jacobian of its argument, and whose execution implements reverse mode AD. In order to do this, ← -J must reflectively access, at runtime, the program in which it is contained, and it must also be possible to apply it to itself. While this offers the possibility of computing higher-order derivatives (in the sense of derivative of the derivative, which we do not consider), it lacks a type-theoretic treatment:

the combinator ← -J is defined in an untyped language. Although [START_REF] Barak | Reverse-mode AD in a Functional Framework: Lambda the Ultimate Backpropagator[END_REF] do mention, for firstorder code, a transformation essentially identical to our ← -D (using so-called backpropagators), the observations that backpropagators may be typed linearly and that ← -D may be directly lifted to higher-order code (as we do in Tab. 3) are original to our work. Finally, let us mention that in this case as well the correctness of ← -J is not illustrated by a mathematical proof of soundness, but by testing an implementation of it (Stalin∇).

At a more theoretical level, [START_REF] Elliott | The simple essence of automatic differentiation[END_REF]] gives a Haskell implementation of backpropagation extracted from a functor D over cartesian categories, with the benefit of disclosing the compositional nature of the algorithm. However, Elliot's approach is still restricted to first-order programs (i.e., computational graphs): as far as we understand, the functor D is cartesian but not cartesian closed, so the higher-order primitives (λ-abstraction and application) lack a satisfactory treatment. This is implicit in Sect. 4.4 of [START_REF] Elliott | The simple essence of automatic differentiation[END_REF], where the author states that he only uses biproduct categories: it is well-known that non-trivial cartesian closed biproduct categories do not exist. 2We already mentioned TensorFlow and PyTorch. It is difficult at present to make a fair comparison with such large-scale differentiable programming frameworks since we are still focused on the conceptual level. Nevertheless, in the perspective of a future implementation, our work is interesting because it would offer a way of combining the benefits of hitherto diverging approaches: the ability to generate modular, optimizable code (TensorFlow) with the possibility of using an expressive language for dynamically-generated computational graphs (PyTorch).

Contents of the paper. Sect. 2 gives a (very subjective) introduction to the notions of AD used in this paper. Sect. 3 informally introduces our main contributions, which will then be detailed in the more technical Sect. 5. Sect. 4 formally specifies the programming language we work with, introducing a simply-typed λ-calculus and a rewriting reduction (Tab. 2) necessary to define and to evaluate our backpropagation transformation ← -D . Sect. 6 applies ← -D to a couple of examples, in particular to a (very simple) recurrent neural network. Sect. 7 concludes with some perspectives.

A CRASH COURSE IN AUTOMATIC DIFFERENTIATION

What follows is an introduction to automatic differentiation for those who are not familiar with the topic. It is extremely partial (the field is too vast to be summarized here) and heavily biased not just towards programming languages but towards the very subject of our work. It will hopefully facilitate the reader in understanding the context and achievements of the paper.

What is automatic differentiation?

Automatic differentiation (or AD) is the science of efficiently computing the derivative of (a restricted class of)

programs [START_REF] Güneş Baydin | Automatic Differentiation in Machine Learning: a Survey[END_REF]. Such programs may be represented as directed acyclic hypergraphs, called computational graphs, whose nodes are variables of type R (the set of real numbers) and whose hyperedges are labelled by functions drawn from some finite set of interest (for example, in the context of neural networks, sum, product and some activation function), with the restriction that hyperedges have exactly one target node and that each node is the target of at most one hyperedge. The basic idea is that nodes that are not target of any hyperedge represent input variables, nodes which are not source of any hyperedge represent outputs, and a hyperedge

x 1 , . . . , x k f -→ y
represents an assignment y := f (x 1 , . . . , x k ), so that a computational graph with n inputs and m outputs represents a function of type R n → R m . An example is depicted in Fig. 1; it represents the function (x 1 , x 2 ) → sin((x 1 -x 2 ) 2 ).

In terms of programming languages, we may define computational graphs as generated by

F , G ::= x | f (x 1 ,. . .,x k ) | let x = G in F | (F ,G)
where x ranges over ground variables and f over a set of real function symbols. The let binders are necessary to represent sharing, a fundamental feature of hypergraphs: e.g., in Fig. 1, node z 1 is shared. The set of real function symbols consists of one nullary symbol for every real number plus finitely many non-nullary symbols for actual functions.

We may write f (G 1 ,. . .,G n ) as syntactic sugar for let x 1 = G 1 in. . . let x n = G n in f (x 1 ,. . .,x n ). Within this introduction (Sect. 2 and 3) we adopt the same notation for a function and the symbol associated with it in the language, so for example r may refer to both a real number and its associated numeral. Also, we use an OCaml-like syntax in the hope that it will help the unacquainted reader parse the examples. From Sect. 4 onward we will adopt a more succinct syntax.

Typing is as expected: types are of the form R k and every computational graph in context x 1 :R,. . .,x n :R ⊢ G : R m denotes a function G : R n → R m (this will be defined formally in Sect. 5.1). Restricting for simplicity to the one-output case, we may say that, as far as we are concerned, the central question of AD is computing the gradient of G at any point r ∈ R n , as efficiently as possible. We remind that the gradient of a differentiable function f : R n → R is defined to be the function ∇f : R n → R n such that, for all r ∈ R,

∇f (r) = (∂ 1 f (r), . . . , ∂ n f (r)) , (1)
where by ∂ i f we denote the partial derivative of f with respect to its i-th argument. Of course, the above question makes sense only if G is differentiable, which is the case if every function symbol represents a differentiable function.

In practice this is not always guaranteed (notoriously, modern neural networks use activation functions which are not differentiable) but this is not actually an issue, as it will be explained momentarily.

Before delving into the main AD algorithms, let us pause a moment on the question of evaluating a computational graph. In terms of hypergraphs, we are given a computational graph G with input nodes x 1 , . . . , x n together with an assignment x i := r i with r i ∈ R for all 1 ≤ i ≤ n. The value G (r 1 , . . . , r n ) is found by progressively computing the assignments w := f (s 1 , . . . , s m ) for each hyperedge z 1 , . . . , z m f -→ w such that the values of all z i are already known.

This is known as forward evaluation and has cost |G | (the number of nodes of G).

In terms of programming languages, forward evaluation corresponds to a standard call-by-value strategy, with values being defined as tuples of numbers. For instance, for G in Fig.

1

let x 1 = 5 in let x 2 = 2 in G - → * let x 1 = 5 in let x 2 = 2 in let x 1 = 3 in let z 2 = z 1 • z 1 in sin z 2 - → * let x 1 = 5 in let x 2 = 2 in let x 1 = 3 in let z 2 = 9 in sin z 2 - → * let x 1 = 5 in let x 2 = 2 in let x 1 = 3 in let z 2 = 9 in 0.412
The operational semantics realizing the above computation will be introduced later. For now, let us mention that the value of a closed computational graph G may be found in O(|G |) reduction steps, where |G | is the size of G as a term, consistently with the cost of forward evaluation.

Forward mode AD

The simplest AD algorithm is known as forward mode differentiation. Suppose that we are given a computational graph (in the sense of a hypergraph) G with input nodes x 1 , . . . , x n and one output node y, and suppose that we want to compute its j-th partial derivative in r = (r 1 , . . . , r n ) ∈ R n . The algorithm maintains a memory consisting of a set of assignments of the form x := (s, t), where x is a node of G and s, t ∈ R (know as primal and tangent), and proceeds as follows:

• we initialize the memory with x i := (r i , 0) for all 0 ≤ i ≤ n, i j, and x j := (r j , 1).

• At each step, we look for a node z 1 , . . . , z k f -→ w such that z i = (s i , t i ) is in memory for all 1 ≤ i ≤ k (there is at least one by acyclicity) and w is unknown, and we add to memory

w := f (s) , k i=1 ∂ i f (s) • t i (2)
where we used the abbreviation s := s 1 , . . . , s m .

This procedure terminates in a number of steps equal to |G | and one may show, using the chain rule for derivatives (which we will recall in a moment), that at the end the memory will contain the assignment y := ( G (r), ∂ j G (r)).

Since the arity k of function symbols is bounded, the cost of computing one partial derivative is O(|G |). Computing the gradient requires computing all n partial derivatives, giving of a total cost of O(n |G |), which is not very efficient since n may be huge (typically, it is the number of weights of a deep neural network, which may well be in the millions).

For example, if G is the computational graph of Fig. 1 and we start with x 1 := (5, 1), x 2 := (2, 0), we obtain z 1 :=

(x 1 -x 2 , 1•1-1•0) = (3, 1), then z 2 := (z 1 •z 1 , z 1 •1+z 1 •1) = (9, 6
) and finally y := (sin(z 2 ), cos(z 2 )•6) = (0.412, -5.467),

which is what we expect since ∂ 1 G (x 1 , x 2 ) = cos((x 1 -x 2 ) 2 ) • 2(x 1 -x 2 ).

Symbolic AD

The AD algorithm presented above is purely numerical, but there is also a symbolic approach. The basic idea of (forward mode) symbolic AD is to generate, starting from a computational graph G with n input nodes and 1 output node, a computational graph -→ D (G) with 2n input nodes and 2 output nodes such that forward evaluation of -→ D (G) corresponds to executing forward mode AD on G, i.e., for all r = r 1 , . . . , r n ∈ R, the output of

- → D (G)(r 1 , 0, . . . , r j , 1, . . . , r n , 0) is ( G (r), ∂ j G (r)).
From the programming languages standpoint, symbolic AD is interesting because:

(1) it allows one to perform optimizations on -→ D (G) which would be inaccessible when simply running the algorithm on G (a typical benefit of frameworks like TensorFlow);

(2) it opens the way to compositionality, with the advantages mentioned in the introduction;

(3) being a (compositional) program transformation rather than an algorithm, it offers a viewpoint from which AD may possibly be extended beyond computational graphs.

Recently, [START_REF] Elliott | The simple essence of automatic differentiation[END_REF]] pointed out how symbolic forward mode AD may be understood in terms of compositionality, thus intrinsically addressing point (2) above. The very same fundamental remark is implicitly used also by [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]. Recall that the derivative of a differentiable function f : R → R is the (continuous) function f ′ : R → R such that, for all r ∈ R, the map a → f ′ (r ) • a is the best linear approximation of f around r . Now, differentiable (resp. continuous) functions are closed under composition, so one may wonder whether the operation (-) ′ is compositional, i.e., whether (д • f ) ′ = д ′ • f ′ . This is notoriously not the case: the so-called chain rule states that

∀r ∈ R, (д • f ) ′ (r ) = д ′ (f (r )) • f ′ (r ). (3)
Nevertheless, there is a slightly more complex construction from which the derivative may be recovered and which has the good taste of being compositional. Namely, define, for f as above,

Df : R × R -→ R × R (r , a) → (f (r ), f ′ (r ) • a)).
We obviously have π 2 Df (x, 1) = f ′ (x) for all x ∈ R (where π 2 is projection on the second component) and we invite the reader to check, using the chain rule itself, that D(д

• f ) = Dд • Df .
The similarity with forward mode AD is not accidental. Indeed, as observed by [START_REF] Elliott | The simple essence of automatic differentiation[END_REF]], forward mode symbolic AD may be understood as a compositional implementation of partial derivatives, along the lines illustrated above.

Formally, we may consider two term calculi for computational graphs, let us call them C and C ′ , defined just as above but based on two different sets of function symbols F ⊆ F ′ , respectively, with F ′ containing at least sum and product, and equipped with partial functions

∂ i : F -→ F ′
for each positive integer i such that, for all f ∈ F of arity k and for all 1 ≤ i ≤ k, ∂ i f is defined and its arity is equal to k. Then, we define a program transformation

- → D from C to C ′ which, on types, is given by - → D (R) := R × R - → D (A × B) := - → D (A) × - → D (B),
and, on computational graphs,

- → D(x) := x - → D(let x = G in F ) := let x = - → D(G) in - → D(F ) - → D((F ,G)) := ( - → D(F ), - → D(G)) - → D (f (x)) := let x = (z, a) in (f (z), k i =1 ∂ i f (z) • a i )
In the last case, f is of arity k and let x = (z, a) in . . . stands for let x 1 = (z 1 ,a 1 ) in . . . let x k = (z k ,a k ) in . . ..

Notice how the case

- → D (f (x))
is the definition of the operator D mutatis mutandis, considering that now the arity k is arbitrary. More importantly, we invite the reader to compare it to the assignment (2) in the description of the forward mode AD algorithm: they are essentially identical. The following is therefore not so surprising:

Proposition 2.1. Suppose that every f ∈ F of arity k corresponds to a differentiable function f : R k → R and that, for all 1 ≤ i ≤ k, ∂ i f is the symbol corresponding to its partial derivative with respect to the i-th input. Then, for every computational graph x : R ⊢ G : R with n inputs, for all 1 ≤ j ≤ n and for all r = r 1 , . . . , r n ∈ R, we have the call-by-value evaluation

let x 1 = (r 1 ,0) in. . . let x j = (r j ,1) in. . . let x n = (r n ,0) in - → D (G) - → * ( G (r),∂ j G (r))
So we obtained what we wanted: evaluating Symbolic AD also provides us with an alternative analysis of the complexity of the forward mode AD algorithm.

- → D (G) in call-
Recall that the length of call-by-value evaluation of computational graphs is linear in the size, so the computation of Proposition 2.1 takes O(

- → D (G) ) steps. Now, by inspecting the definition of - → D, it is immediate to see that every syntactic
construct of G is mapped to a term of bounded size (remember that the arity k is bounded). Therefore,

- → D (G) = O(|G |),
which is exactly the cost of computing one partial derivative in forward mode. In other words, the cost becomes simply the size of the output of the program transformation.

Notice how Proposition 2.1 rests on a differentiability hypothesis. As mentioned above, this is not truly fundamental.

Indeed, observe that -→ D (G) is always defined, independently of whether the function symbols it uses are differentiable. This is because, in principle, the maps ∂ i : F → F ′ are arbitrary and the symbol ∂ i f may have nothing to do with the i-th partial derivative of the function that the symbol f represents! Less unreasonably, we may consider "mildly" non-differentiable symbols and associate with them "approximate" derivatives: for example, the rectified linear unit ReLU(x) defined as if x<0 then 0 else x may be mapped by ∂ 1 to Step(x) defined as if x<0 then 0 else 1, even though technically the latter is not its derivative because the former is not differentiable in 0. Formally, Step is called a subderivative of ReLU. Proposition 2.1 easily extends to subderivatives and, therefore, AD works just as well on non-differentiable functions, computing subgradients instead of gradients, which is perfectly acceptable in practice.

This remark applies also to backpropagation and to all the results of our paper, although for simplicity we prefer to stick to the more conventional notion of gradient, and thus work under a differentiability hypothesis.

Reverse mode AD, or backpropagation

A more efficient way of computing gradients in the many inputs/one output case is provided by reverse mode automatic differentiation, from which the backpropagation (often abbreviated as backprop) algorithm derives. As usual, we are 3 Technically, one may consider the calculi C and C ′ as cartesian 2-multicategories: types are objects, computational graphs with n inputs and one output are n-ary multimorphisms and evaluation paths are 2-arrows. Then, -→ D is readily seen to be a morphism of cartesian 2-multicategories. We will not develop such categorical considerations in this paper, we will content ourselves with mentioning them in footnotes.

given a computational graph (seen as a hypergraph) G with input nodes x 1 , . . . , x n and output node y, as well as r = r 1 , . . . , r n ∈ R, which is the point where we wish to compute the gradient. The backprop algorithm maintains a memory consisting of assignments of the form x := (r , α), where x is a node of G and r , α ∈ R (the primal and the adjoint), plus a boolean flag with values "marked/unmarked" for each hyperedge of G, and proceeds thus: initialization: the memory is initialized with x i := (r i , 0) for all 1 ≤ i ≤ n, and the forward phase starts;

forward phase: at each step, a new assignment z := (s, 0) is produced, with s being computed exactly as during forward evaluation, ignoring the second component of pairs in memory (i.e., s is the value of node z); once every node of G has a corresponding assignment in memory, the assignment y := (t, 0) is updated to y := (t, 1), every hyperedge is flagged as unmarked and the backward phase begins (we actually know that t = G (r), but this is unimportant); backward phase: at each step, we look for an unmarked hyperedge z 1 , . . . , z k f -→ w such that all hyperedges having w among their sources are marked. If no such hyperedge exists, we terminate. Otherwise, assuming that the memory contains w := (u, α) and z i := (s i , β i ) for all 1 ≤ i ≤ k, we update the memory with the assignments

z i := (s i , β i + ∂ i f (s) • α) (where s := s 1 , . . . , s k ) for all 1 ≤ i ≤ k
and flag the hyperedge as marked.

One may prove that, when the backprop algorithm terminates, the memory contains assignments of the form

x i := (r i , ∂ i G (r)) for all 1 ≤ i ≤ n, i.e.
, the value of each partial derivative of G in r is computed at the corresponding input node, and thus the gradient may be obtained by just collecting such values. Let us test it on an example. Let G be the computational graph of Fig. 1 and let r 1 = 5, r 2 = 2. The forward phase terminates with x 1 := (5, 0), x 2 := (2, 0), z 1 := (3, 0), z 2 := (9, 0), y := (0.412, 1). From here, the backward phase updates z 2 := (9, cos(z 2 ) • 1) = (9, -0.911), then

z 1 := (3, z 1 • -0.911 + z 1 • -0.911) = (3, -5.467) and finally x 1 := (5, 1 • -5.467) = (5, -5.467) and x 2 := (2, -1 • -5.467) = (2, 5.467), as expected since ∂ 2 G = -∂ 1 G .
Compared to forward mode AD, the backprop algorithm may seem a bit contrived (it is certainly harder to understand why it works) but the gain in terms of complexity is considerable: the forward phase is just a forward evaluation; and, by construction, the backward phase scans each hyperedge exactly once performing each time a constant number of operations. So both phases are linear in |G |, giving a total cost of O(|G |), like forward mode. Except that, unlike forward mode, a single evaluation now already gives us the whole gradient, independently of the number of inputs!

Symbolic backpropagation and the compositionality issue

The symbolic methodology we saw for forward mode AD may be applied to the reverse mode too: given a computational graph G with n inputs and one output, one may produce a computational graph bp(G) with n + 1 inputs and 1 +n outputs such that, for all r = r 1 , . . . ,

r n ∈ R, the forward evaluation of bp(G)(r, 1) has output ( G (r), ∇ G (r)). Moreover, bp(G) = O(|G |).
A formal definition of the bp transformation will be given in Sect. 5.1. Let us look at an example, shown in Fig. 2.

First of all, observe that bp x 1 ,x 2 ,a (G) contains a copy of G, marked in blue. This corresponds to the forward phase. The nodes marked in red correspond to the backward phase. Indeed, we invite to reader to check that the forward evaluation of bp x 1 ,x 2 ,a (G) with x 1 := 5, x 2 := 2 and a := 1 matches exactly the steps of the backprop algorithm as exemplified in Sect. 2.4, with node b (resp. c, x ′ 1 , x ′ 2 ) corresponding to the second component of the value of node z 2 (resp. z 1 , x 1 , x 2 ). (The nodes v, c ′ and c ′′ are just intermediate steps in the computation of b and c which are implicit in the numerical description of the algorithm and are hidden in syntactic sugar).

x 1 x 2 z 1 z 2 y - • sin a b v c ′ c ′′ c x ′ 1 x ′ 2 cos • • • + 1 • -1 • let z 1 =x 1 -x 2 in let z 2 =z 1 •z 1 in let b=(cos z 2 )•a in let c=z 1 •b+z 1 •b in (sin z 2 ,(1•c,-1•c))
Fig. 2. The computational graph bp x 1 , x 2 , a (G) where G is in Fig. 1, and its corresponding term.

Rather than examining the details of the definition of bp, let us observe at once that, from the standpoint of programming languages, it suffers from a serious defect: unlike -→ D, it is not compositional. Indeed, in order to define bp(let x = G in F ), we need to know and exploit the inner structure of bp(F ) and bp(G), whereas from the definition of -→ D given above it is clear that no such knowledge is required in that case, i.e., -→ D (F ) and -→ D (G) are treated as "black boxes". Our way of understanding previous work on the subject [START_REF] Elliott | The simple essence of automatic differentiation[END_REF][START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]] is that it was all about making symbolic backprop compositional. This is the topic of the next Section.

OUR APPROACH TO COMPOSITIONAL BACKPROPAGATION

As mentioned above, the key to modular and efficient differentiable programming is making symbolic backprop (the bp transformation) compositional. We will show that this may be achieved by considering a programming language with a notion of linear negation. The goal of this section is to explain why negation and why linear.

Let us start with by looking at an extremely simple example, which is just the composition of two unary functions:

G := let z = f x in д z (4)
As a hypergraph, G has three nodes x, y, z, of which x is the input and y the output (corresponding to the root of G)

and two edges x f -→ z and z д -→ y. Note that, since G has only one input, its gradient is equal to its derivative and forward and reverse mode AD have the same complexity. This is why the example is interesting: it shows the difference between the two algorithms by putting them in a context where they must perform essentially the same operations. In the sequel, we set h := G and we denote by f ′ , д ′ and h ′ the derivatives of f , д and h, respectively.

For what concerns forward mode, we invite the reader to check that

- → D (G) = let z = let (v,a)= x in (f v, (f ′ v) • a) in let (w, b) = z in (д w, (д ′ w) • b)
We may simplify this a bit by observing that (renaming w as z)

- → D (G){(x,a)/x} - → * let z = f x in let b = (f ′ x)•a in (д z, (д ′ z)•b)
On the other hand, applying the definition of Sect. 5.1, we get

bp x,a (G) = let z = f x in let b = (д ′ z)•a in (д z,(f ′ x)•b)
Note that, if we substitute r ∈ R for x and 1 for a, in both cases we obtain (h(r ), h ′ (r )) in the same number of steps, as expected. However, the order in which the derivatives are computed relative to the order of the composition д • f is different: it is covariant in forward mode, contravariant in reverse mode. This corresponds precisely to the behavior of the two algorithms:

• in forward mode, we start with x := (r , 1), from which we infer z := (f (r ), f ′ (r )), from which we infer

y := (д(f (r )), д ′ (f (r )) • f ′ (r ));
• in reverse mode, the forward phase leaves us with x := (r, 0), z := (f (r ), 0), y := (д(f (r )), 1), at which point the backward phase proceeds back from the output towards the input, inferring first z := (f (r ), д ′ (f (r ))) and then

x := (r , f ′ (r ) • д ′ (f (r ))).
A lesson we learn from this example, in the perspective of compositionality, is that both algorithms may be seen as mapping the subprograms f and д, which have one input and one output, to two subprograms

- → f and - → д (or ← - f and ← - д )
having two inputs and two outputs, but then assembling them rather differently:

forward mode reverse mode

- → f - → д r h(r ) 1 h ′ (r ) ← - f ← - д r h(r ) h ′ (r ) 1
The picture on the right suggested to [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]] the idea of solving the compositionality issue via continuations:

drawing inspiration from "There and Back Again" [START_REF] Danvy | There and Back Again[END_REF], the blocks ←f and ←д are seen as function calls in the CPS transform of G, so that the forward phase takes place along the call path, while the backward phase takes place along the return path. However, in order for their approach to work, [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF] must use references,

i.e., the memory maintained by the backprop algorithm is explicitly present and is manipulated as described in Sect. 2.4.

Moreover, since the memory is updated after the return from each function call, they must actually resort to delimited continuations. On the whole, although they do succeed in presenting reverse mode AD as a compositional program transformation, the work of [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]] is closer to an (elegant!) implementation of the backprop algorithm in a functional language with references than to an abstract, purely compositional description of its dynamics.

Let us focus, instead, on the idea of contravariance. The archetypal contravariant operation is negation. For a (real)

vector space A, negation corresponds to the dual space A ⊸ R, which may be generalized to A ⊥ E := A ⊸ E for an arbitrary space E, although in fact we will always take E = R d for some d ∈ N. For brevity, let us keep E implicit and simply write A ⊥ . There is a canonical way, resembling CPS, of transforming a differentiable function

f : R → R with derivative f ′ into a function D r f : R × R ⊥ → R × R ⊥ from
which the derivative of f may be extracted. Namely, let, for

all x ∈ R and x * ∈ R ⊥ , D r f (x, x * ) := f (x), λa.x * (f ′ (x) • a) , (5)
where we are using λ-notation with the standard meaning. We let the reader verify that, if we suppose

E = R, then for all x ∈ R, (π 2 D r f (x, I ))1 = f ′ (x)
, where π 2 is projection of the second component and I : R → R is the identity function (which is obviously linear). More importantly, D r is compositional: for all x ∈ R and x * ∈ R ⊸ R, we have 4

D r д(D r f (x, x * )) = D r д(f (x), λa.x * (f ′ (x) • a)) = (д(f (x)), λb.(λa.x * (f ′ (x) • a))(д ′ (f (x)) • b)) = (д(f (x)), λb.x * (f ′ (x) • (д ′ (f (x)) • b)) = (д(f (x)), λb.x * ((д ′ (f (x)) • f ′ (x)) • b)) = ((д • f )(x), λb.x * ((д • f ) ′ (x) • b)) = D r (д • f )(x, x * ).
This observation may be generalized to maps f : R n → R: for all x ∈ R n and

x * = x * 1 , . . . , x * n ∈ R ⊥ , ← - D(f )(x; x * ) := f (x), λa. n i=1 x * i (∂ i f (x) • a) .
In the AD literature, the x * i are called backpropagators [Pearlmutter and [START_REF] Barak | Reverse-mode AD in a Functional Framework: Lambda the Ultimate Backpropagator[END_REF]. Obviously

← - D(f ) : (R × R ⊥ ) n → R × R ⊥
and we invite the reader to check that, if we take E = R n , we have

(π 2 ← - D(f )(x; ι 1 , . . . , ι n ))1 = ∇f (x),
where, for all 1 ≤ i ≤ n, ι i : R → R n is the injection into the i-th component, i.e., ι i (x) = (0, . . . , x, . . . , 0) with zeros everywhere except at position i, which is a linear function. Moreover, ← -D is compositional. 5 This leads to the definition of a compositional program transformation

← - D (Tab. 3) which verifies ← - D(R) = R × R ⊥
and which, applied to our example (4), gives

← - D(G) = let z = let (v,v * ) = x in (f v, fun b -> v * ((f ′ v) • b)) in let (w, w * ) = z in (д w, fun a -> w * (д ′ (w) • a))
where w * , v * : R ⊥ so that both fun b -> v * ((f ′ v) • b)) and fun a -> w * ((д ′ w) • a) have also type R ⊥ . Notice the resemblance of ← -D(G) and -→ D(G): this is not an accident, both are defined in a purely compositional way (in particular, ← -

D(let x = H in F ) = let x = ← - D(H ) in ← - D(F ))
, abiding by the "black box" principle) and the only non trivial case is when they are applied to function symbols. Moreover, we have

← - D(G){(x,x * )/x} - → * let z = f x in (д z, fun a -> let b = (д ′ z)•a in x * ((f ′ x)•b))
which, albeit of different type, is essentially identical to bp x,a (G). More precisely, if we write

let b = (д ′ z)•a in (f ′ x)•b as F , then we have bp x,a (G) =let z=f x in (д z,F ) and ← - D(G) {(x,x * )/x} - → * let z=f x in (д z, fun a -> x * F )
Let us now explain why negation must be linear. The above example is too simple for illustrating this point, so from now on let G be the computational graph of Fig. 1. Applying the definition of Tab. 3 and simplifying, we obtain that ← -D (G) is equal to: 4 The equation in the second line uses both commutativity and associativity of product, but only the latter is really necessary: by replacing f ′ (x ) • a with a • f ′ (x ) in Eq. 5 one can check that commutativity is not needed. The former notation reflects that this is actually a linear application: in general, if f : A → B, then f ′ (x ) : A ⊸ B and a : A. When A = B = k with k a ring (commutative or not), k ⊸ k k and linear application becomes the product of k, so the notation a • f ′ (x ) makes sense and backprop has in fact been applied to non-commutative rings [START_REF] Isokawa | Quaternion Neural Network and Its Application[END_REF][START_REF] Pearson | Back Propagation in a Clifford Algebra[END_REF]]. In the general case, it makes no sense to swap function and argument and it is unclear how backprop would extend. 5 Technically, functions of type R n → R for varying n form what is known as a cartesian operad, or clone, and ← -D is a morphism of such structures. A bit more explicitly, one may compose f : R n → R with д : R m+1 → R by "plugging" f into the i-th coordinate of д, forming д

• i f : R n+m → R, for any 1 ≤ i ≤ m + 1; the operation ← - D preserves such compositions. let z 2 = let z 1 = let (v 2 , v * 2 ) = x 2 in let (v 1 , v * 1 ) = x 1 in (v 1 -v 2 , fun c -> v * 1 (1•c)+v * 2 (-1•c)) in let (w 1 ,w * 1 ) = z 1 in (w 1 •w 1 , fun b -> w * 1 (w 1 •b)+w * 1 (w 1 •b)) in let (w 2 , w * 2 ) = z 2 in (sin w 2 , fun a -> w * 2 ((cos w 2 )•a)) - → * let (v 2 , v * 2 ) = x 2 in let (v 1 , v * 1 ) = x 1 in let (z 1 ,z * 1 ) = (v 1 -v 2 , fun c -> v * 1 (1•c)+v * 2 (-1•c)) in let (z 2 , z * 2 ) = (z 1 •z 1 , fun b -> z * 1 (z 1 •b)+z * 1 (z 1 •b)) in (sin z 2 , fun a -> z * 2 ((cos z 2 )•a))
There is a potential issue here, due to the presence of two occurrences of z * 1 (highlighted in brown) which are matched against the function corresponding to the derivative of x 1 -x 2 (also highlighted in brown, let us denote it by F ). Notice that such a derivative is present only once in bp x 1 ,x 2 ,a (G): it corresponds to the rightmost nodes of Fig. 2 (more precisely, the abstracted variable c corresponds to node c itself, whereas v * 1 and v * 2 correspond to nodes x ′ 1 and x ′ 2 , respectively). Therefore, duplicating F would be a mistake in terms of efficiency.

The key observation here is that z * 1 is of type R ⊥ , i.e., it is a linear function (and indeed, F is linear: by distributivity of product over sum, c morally appears only once in its body). This means that, for all t, u : R, we have z * 1 t +z * 1 u = z * 1 (t +u). In the λ-calculus we consider (Sect. 4), this becomes an evaluation step oriented from left to right, called linear factoring still computes ∇ t with the same efficiency as the evaluation of t! Indeed, the definition of

← - D immediately gives us that if t - → * u in p steps, then ← - D (t) - → * ← - D (u) in O(p)
steps (point 2 of Lemma 5.4). 6 But since t has ground type and ground free variables, eliminating all higher-order redexes gives t -→ * G for some computational graph G, hence ← -

D (t)
evaluates to (essentially) bp(G). So ← -D (t) computes the gradient of t = G (remember that the semantics is invariant under evaluation) with a cost equal to O(|G |) plus the cost of the evaluation t -→ * G, which is the best we can expect in general.

To conclude, we should add that in the technical development we actually use Accattoli's linear substitution calculus

[Accattoli 2012],
which is a bit more sophisticated than the standard simply-typed λ-calculus. This is due to the presence of linear negation, but it is also motivated by efficiency, which is the whole point of backpropagation. To be taken seriously, a proposal of using functional programming as the foundation of (generalized) AD ought to come with a thorough complexity analysis ensuring that we are not losing efficiency in moving from computational graphs to λ-terms. Thanks to its tight connection with abstract machines and reasonable cost models [START_REF] Accattoli | Distilling Abstract Machines[END_REF]],

6 Morally, the simply-typed λ-calculus augmented with a set F of function symbols is the free cartesian semi-closed 2-multicategory on F (semi-closed in the sense of [START_REF] Hyland | Classical lambda calculus in modern dress[END_REF]). Therefore, once ← -D is defined on F, it automatically extends to a morphism of such structures. In particular, it functorially maps evaluations (which are 2-arrows) to evaluations.

ultimately owed to its relationship with Girard's proof nets [START_REF] Accattoli | Proof Nets and the Linear Substitution Calculus[END_REF]] (a graphical representation of linear logic proofs which may be seen as a higher order version of computational graphs), the linear substitution calculus is an ideal compromise between the abstractness of the standard λ-calculus and the concreteness of implementations, and provides a solid basis for such an analysis.

THE LINEAR SUBSTITUTION CALCULUS

Terms and Types. Since the linear factoring rule ( 18) is type-sensitive (as mentioned above, it is unsound in general), it is convenient to adopt a Church-style presentation of our calculus, i.e., with explicit type annotations on variables.

The set of types is generated by the following grammar:

A, B, C ::= R | A × B | A → B | R ⊥ d (simple types)
where R is the ground type of real numbers. The negation R ⊥ d corresponds to the linear implication (in the sense of linear logic [START_REF] Girard | Linear Logic[END_REF]]) R ⊸ R d . However, in order to keep the calculus as simple as possible, we avoid introducing a fully-fledged bilinear application (as for example in the bang-calculus [START_REF] Ehrhard | The Bang Calculus: an untyped lambda-calculus generalizing call-by-name and call-by-value[END_REF]) and opt instead for just a negation operator and dedicated typing rules. We may omit the subscript d in R ⊥ d if clear from the context or unimportant.

An annotated variable is either x !A (called exponential variable) with A any type, or x R (called linear variable): the writing x (!)A stands for one of the two annotations (in the linear case A = R). The grammar of values and terms is given by mutual induction as follows, with x (!)A varying over the set of annotated variables, r over the set of real numbers R and f over a finite set F of function symbols over real numbers containing at least multiplication (noted in infix form t • u):

v ::= x (!)A | r | λx (!)A .t | ⟨v 1 , v 2 ⟩ (values) (6) t, u ::= v | tu | ⟨t, u⟩ | t[ x !A , y !B := u] | t[x (!)A := u] | t + u | f (t 1 , . . . , t k ) (terms) (7)
Since binders contain type annotations, bound variables will never be annotated in the sequel. In fact, we will entirely omit type annotations if irrelevant or clear from the context. Terms of the form r are called numerals. The term t[x := u]

(and its binary version t[⟨x, y⟩ := u]) may be understood as the more familiar let x = u in t used in the previous informal sections.

We denote by Λ ⊥ (F ) the set of terms generated by the above grammar. We consider also the subset Λ(F ) of terms obtained by discarding linear negation and linear variables.

We denote by |t | the size of a term t, i.e., the number of symbols appearing in t. We denote by fv(t) the set of free variables of t, abstractions and explicit substitutions being the binding operators. As usual, α-equivalent terms are treated as equal. A term t is closed if fv(t) = ∅. In the following, we will use boldface metavariables to denote sequences:

x will stand for a sequence of variables x 1 , . . . , x n , t for a sequence of terms t 1 , . . . , t n , etc. The length of the sequences will be left implicit, because either clear from the context or irrelevant.

We use n-ary products ⟨t 1 , . . . , t n-1 , t n ⟩ as syntactic sugar for ⟨t 1 , . . . ⟨t n-1 , t n ⟩ . . . ⟩, and we define Euclidean types by R d := R × (. . . R × R . . .). It will be useful to denote a bunch of sums without specifying the way these sums are associated. The notation

Γ ⊢ (z) t • u : R Γ ⊢ z 0 : R d Γ ⊢ (z) t : R d Γ ⊢ (z) u : R d Γ ⊢ (z) t + u : R d
Table 1. The typing rules. In the pairing and sum rules, either all three sequents have z, or none does.

will denote such a bunch for I a finite set. In case I is a singleton, the sum denotes its unique element. An empty sum of type R d stands for 0, . . . , 0 , which we denote by 0. Of course this notation would denote a unique term modulo associativity and commutativity of + and neutrality of 0, but we do not need to introduce these equations in the calculus.

The typing rules are in Tab. 1. The meta-variables Γ, ∆ vary over the set of typing contexts, which are finite sequences of exponential type annotated variables without repetitions. There are two kind of sequents: Γ ⊢ t : A and Γ ⊢ z t : R d .

In this latter d ∈ N and z is linear type annotated variable which occurs free linearly in t. The typing rules define what "occurring linearly" means, following the standard rules of linear logic. 7 The writing Γ ⊢ (z) t : A stands for either Γ ⊢ t : A or Γ ⊢ z t : A, and in the latter case A = R d for some d.

Contexts. We consider one-hole contexts, or simply contexts, which are defined by the above grammar (6) restricted to the terms having exactly one occurrence of a specific variable {•}, called the hole. Meta-variables C, D will range over the set of one-hole contexts. Given a context C and a term t we denote by C{t } the substitution of the hole in C by t allowing the possible capture of free variables of t. A particular class of contexts are the substitution contexts which have the form of a hole followed by a stack (possibly empty) of explicit substitutions:

{•}[p 1 := t 1 ] . . . [p n := t n ]
with each p i a variable or a pair of variables. Meta-variables α, β will range over substitution contexts. If α is a substitution context, we will use the notation tα instead of α {t }.

Rewriting rules. The reduction relation is given in Tab. 2 divided in three sub-groups: 8), (9), ( 10), ( 11), ( 12), ( 13), ( 14), (15)} evaluation rules, η := {( 16), (17)} extensional rules, ℓ := {(18)} linear factoring.

β := {(
In case one wants to consider numeric computations (other than sum and products), then of course one must also include suitable reduction rules associated with the function symbols:

f (r 1 α 1 , . . . , r n α n ) - → f (r 1 , . . . , r n )α 1 . . . α n (24) (λx .t)αu - → t[x := u]α (8) s[⟨x, y⟩ := ⟨t, u⟩ α] - → s[x := t][y := u]α (9) C{x }[x !A := vα] - → C{v}[x !A := v]α (10) t[x !A := vα] - → tα if x fv(t) (11) t[x R := vα] - → t {v/x } α (12) rα + qβ - → r + qα β (13) ⟨t 1 , t 2 ⟩ α + ⟨u 1 , u 2 ⟩ β - → ⟨t 1 + u 1 , t 2 + u 2 ⟩ α β (14) rα • qβ - → rqα β (15)
(a) β -rules. In (10), ( 11) and (12), v is a value. In (10), C is an arbitrary context not binding x . We write (10) n to refer to the instance of (10) in which v is a numeral.

t -→ λy.ty

(16) t[x := u, u ′ ] - → t[x := y, y ′ ][ y, y ′ := u, u ′ ] (17)
(b) η-rules. In (16) t has an arrow type or R ⊥ . The new variables on the right-hand side of both rules are fresh.

(x R ⊥ αt)β + (x R ⊥ α ′ t ′ )β ′ - → x R ⊥ (t + t ′ )α βα ′ β ′ (18)
(c) linear factoring (ℓ-rule for short), where we suppose that none of α, β, α ′ , β ′ binds x .

t

[x := u][y := w] ≡ t[y := w][x := u] if y fv(u) and x fv(w) (19) t[x := u][y := w] ≡ t[x := u[y := w]] if y fv(t) (20) t[x !A := u] ≡ t {y/x } [x !A := u][y !A := u] (21) (s□t)[x := u] ≡ s[x := u]□t if x fv(t) (22) (s□t)[x := u] ≡ s□(t[x := u]) if x fv(s) (23)
(d) Structural equivalence. In (21), t {y/x } denotes t in which some (possibly none) occurrences of x are renamed to a fresh variable y.

In (22), ( 23) the writing s□t stands for either st or s + t or s • t or ⟨s, t ⟩.

Table 2. The reduction and the structural equivalence relations, where we suppose the usual convention that no free variable in one term can be captured in the other term of a relation.

The rule (8) transforms a λ-calculus application into an explicit substitution. The difference between the two is that the latter commutes over terms by the structural equivalence defined in Tab. 2d, while the former does not. Rule (9) deconstructs a pair, while rule (10) implements a "micro-step" substitution, closer to abstract machines [START_REF] Accattoli | Distilling Abstract Machines[END_REF]]. The special case in which v is a numeral is referred to as (10) n . Rule (11) implements garbage collection, and rule (12) linear substitution. The rules ( 16) and ( 17) are standard instances of η-expansion rules. They are useful in the proof of Theorem 5.6. Rule (18) has already been discussed.

Notice that Λ(F ) is a standard linear explicit substitution calculus encompassing both call-by-need and call-byvalue [START_REF] Accattoli | Distilling Abstract Machines[END_REF]. For instance, the usual by-value β-rule (λx .t)v -→ t {v/x } is derivable. In this respect, the reader may think of Λ(F ) as nothing but the plain simply-typed λ-calculus, and consider explicit substitutions as needed merely to represent computational graphs (which are obtained by restricting to the ground type R only). The situation is different in Λ ⊥ (F ), in which linearity plays a key role for expressing backpropagation.

Given any X ⊆ β ∪ η ∪ ℓ, we denote by X --→ the context closure of the union of the reduction relations in X , for any context C:

C{t } X --→ C{u}, whenever t X --→ u.
This means that we do not consider a specific evaluation strategy (call-by-value, call-by-name etc. . . ), in order to be as general as possible and to allow a future analysis concerning a more efficient operational semantics.

A term t is a X -normal form if there is no term u with t Structural equivalence ≡ is the smallest equivalence relation containing the context closure of the rules ( 19)-( 23) in Tab. 2d. Morally, structural equivalence relates terms which would correspond to the same state of an abstract machine implementing the calculus, in which explicit substitutions represent pointers to memory. We refer to [START_REF] Accattoli | Distilling Abstract Machines[END_REF][START_REF] Accattoli | Environments and the complexity of abstract machines[END_REF] for more details. The crucial property of ≡ is that it is a (strong) bisimulation with respect to -→, which means in particular that it may always be postponed to the end of an evaluation (Proposition 4.3).

The following properties are standard and we give them without proof.

Proposition 4.1 (Subject reduction). If t -→ u or t ≡ u and Γ ⊢ (z) t : A, then Γ ⊢ (z) u : A.

Lemma 4.2 (≡ is a strong bisimulation). Let ι be any reduction rule and let

t ′ ≡ t ι - → u, then there exists t ′ ι - → u ′ such that u ′ ≡ u.
Proposition 4.3 (postponement of ≡). Let X be any subset of the reduction rules in Tab. 2 (including the variant (10) n ) and let t ( X --→ ∪ ≡) * u with k X -steps, then there exists u ′ such that t X --→ k u ′ ≡ u.

Proposition 4.4 (values). Given a closed term t of type A, if t is a β-normal form, then it is a value.

Proposition 4.5 (Weak normalization). For every term t, and every set X ⊆ βℓ, there exists a X -normal form u such that t X --→ * u.

The βℓ-rewriting enjoys also strong normalization, even modulo ≡, but the proof is more involved and uninteresting for our purposes, so we omit it. The strong normalization is however immediate if we restrict the contraction rule (10) to numerals, a property which will be useful in the sequel.

Lemma 4.6. For any t (10) n (11)(14)(13)( 15)

----------------→ * u, the number of steps in the sequence is O(|t |).

Denotational semantics. The cartesian closed category of sets and functions gives a denotational model for this calculus. Types are interpreted as sets, as follows:

R := R A → B := set of functions from A to B A × B := A × B R ⊥ d := set of linear maps from R to R d
Notice that the restriction to linear functions in R ⊥ d is such that rule ( 18

* v ′ . Then v = v ′ .
Proof. It is not hard to show that the interpretation is injective on tuples of numerals, i.e., for v, v

′ : R d values, v = v ′ implies v = v ′ .
The statement then follows from Props. 4.4 and 4.7.

□

From now on, we will suppose that:

(⋆) all function symbols in F are associated bywith differentiable maps on real numbers.

As mentioned at the end of Sect. 2.3, this hypothesis is essentially cosmetic, in that it allows us to use actual gradients and to avoid the more technical notion of subgradient (see also Sect. 7).

Any term x !R 1 , . . . , x !R n ⊢ t : R is denoted by an n-ary map t over R. Since differentiable functions compose, then (⋆) implies that t is also differentiable, if t contains only variables of type R. In the general case, one can use Prop. 4.5 in order to β-reduce t into a β-normal form u containing only variables of type R. Then by Prop. 4.7 t = u and so t is differentiable. This justifies the following notation, given x !R 1 , . . . , x !R n ⊢ t : R and a vector r ∈ R n :

∇t(r) := ∂ 1 t (r), . . . , ∂ n t (r) . (25)
Section 5 gives an efficient way of computing ∇t(r) based on the syntactic structure of t.

THE BACKPROPAGATION TRANSFORMATION

Let us fix two sets of function symbols F , F ′ such that F ⊆ F ′ , together with partial functions ∂i : F -→ F ′ , for each positive integer i such that, for all f ∈ F of arity k and for all 1 ≤ i ≤ k, ∂ i f is defined and its arity is equal to k.

In addition to the hypothesis (⋆) in Sect. 4, we also suppose:

(⋆⋆) ∂ i f := ∂ i f For any d ∈ N, Tab. 3 defines a program transformation ← - D d from Λ(F ) to Λ ⊥ (F ′ )
, called the reverse gradient relative to R d . Given a term x !R ⊢ t : R, Cor. 5.7 proves that ← -D d (t) computes the gradient of t in at most O(m + |G |) steps, where m is the cost of evaluating t to a computational graph G. Moreover, Cor. 5.7 is a consequence of Th. 5.6, proving that actually one can reduce ← -D d (t) to a term expressing the backpropagation algorithm applied to any computational graph G β-equivalent to t (indeed, to a single λ-term t one can associate different computational graphs, with different sizes and with different degrees of sharing). Sect. 5.1 sets the framework needed to state and prove our results. Grammar (26) defines the notion of ground term corresponding to a computational graph and Def. 5.1 gives the computational graph associated with the backpropagation applied to ground terms. Prop. 5.3 formally proves the soundness of this algorithm. Sect. 5.2 then moves to the more general case of λ-terms, giving the soundness of our reverse gradient transformation ← -

D d . ← - D d (R) := R × R ⊥ d ← - D d (A → B) := ← - D d (A) → ← - D d (B) ← - D d (A × B) := ← - D d (A) × ← - D d (B) (a)
The action of the transformation on types.

← -

D d (x !A ) := x ! ← - D d (A) ← - D d (λx !A .t) := λx ! ← - D d (A) . ← - D d (t) ← - D d (tu) := ← - D d (t) ← - D d (u) ← - D d (⟨t, u⟩) := ← - D d (t), ← - D d (u) ← - D d (t[ x !A , y !B := u]) := ← - D d (t)[ x ! ← - D d (A) , y ! ← - D d (B) := ← - D d (u)] ← - D d (t[x !A := u]) := ← - D d (t)[x ! ← - D d (A) := ← - D d (u)] ← - D d (r ) := r , λa R .0 ← - D d (t + u) := x + y, λa R .(x * a + y * a) [ x !R , x * !R ⊥ d := ← - D d (t)][ y !R , y * !R ⊥ d := ← - D d (u)] ← - D d (f (t)) := f (x) , λa R . k i=1 x * i (∂ i f (x) • a) [ x !R , x * !R ⊥ d := ← - D d (t)] (b)
The action of the transformation on terms. In the definition of ← -D (f (t)), the sequences t, x, x * have all length k equal to the arity of f and the notation

[ x !R , x * !R ⊥ d := ← - D (t)] stands for [ x !R 1 , x * 1 !R ⊥ d := ← - D (t 1 )] • • • [ x !R k , x * k !R ⊥ d := ← - D (t k )].
As mentioned in Sect. 3, the variables with superscript * correspond to backpropagators in AD terminology [START_REF] Barak | Reverse-mode AD in a Functional Framework: Lambda the Ultimate Backpropagator[END_REF]. may be arbitrary (it plays a role only in the last remark of Sect. 5.2) is a crucial feature allowing its compositionality, in contrast with the definition of bp x,a (G) which has to refer to x containing the free variables in G (see the case bp x,a (f (y)) in Def. 5.1). This being said, we henceforth omit the index d.

Backpropagation on Computational Graphs

We restrict Λ(F ) to terms of type R not containing higher types, deemed ground terms, as follows:

F , G ::= x !R | F [x !R := G] | f (x !R 1 , . . . , x !R k ) | r | F + G (26) A term f (G 1 , . . . , G k ) is considered as syntactic sugar for f (x 1 , . . . , x k )[x 1 := G 1 ] . . . [x k := G k ].
Fig. 1 and2 give examples of ground terms with the associated computational graph. Notice that the type system of Tab. 1 assigns to any ground term G a type judgment x !R 1 , . . . , x !R n ⊢ G : R for a suitable set of variables, so ∇G(r) is defined by ( 25) and the hypothesis (⋆) and (⋆⋆).

We now define the transformation bp implementing symbolic backpropagation, as described on hypergraphs e.g. in [Van Iwaarden 1993, Sect. 3]. We first introduce the following notation, evaluating some trivial sums on the fly: given two ground terms F 0 , F 1 ,

F 0 ⊕ F 1 :=        F i if F i-1 = 0, F 0 + F 1 otherwise Definition 5.1. Let Γ = x !R 1 , . . . , x !R n .
Given a ground term of type Γ ⊢ G : R and a fresh variable a R , we define the term

Γ, a !R ⊢ bp x,a (G) : R × R n
by induction on G, as follows. The definition is based on the inductive invariant that

bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β
where α and β are substitution contexts and a does not appear free in G 0 or in β.

• bp x,a (x i ) := x i , 0, . . . , a, . . . , 0 , where a appears at the i-th position in the tuple.

• bp x,a (r ) := r, 0, . . . , 0 .

• Let f be of arity k and y a subsequence of length k of x. Then, bp x,a (f (y)) := f (y), 0, . . . , 0,

∂ 1 f (y) • a, 0, . . . , 0, ∂ k f (y) • a, 0, . . . , 0 ,
where the non-zero terms in the tuple are at the positions determined by y within x.

• Let G = F ′ [z !R := F ′′ ] and suppose that (with b !R a new variable distinct from a !R ) bp x,z,a (F ′ ) = F ′ 0 , F ′ 1 , . . . , F ′ n , H α ′ β ′ , bp x,b (F ′′ ) = F ′′ 0 , F ′′ 1 , . . . , F ′′ n α ′′ β ′′ .
Notice that for every i ≤ n we can suppose by renaming that the set of variables of F ′ i bound by α ′ , β ′ is disjoint from the set of variables of F ′′ i bound by α ′′ , β ′′ . Also, notice that F ′′ 0 has type x !R ⊢ F ′′ 0 : R. So we can define:

bp x,a (F ′ [z !R := F ′′ ]) :=        F ′ 0 , F ′ 1 , . . . , F ′ n α ′ β ′ [z !R := F ′′ 0 ]β ′′ if H = 0, F ′ 0 , F ′ 1 ⊕ F ′′ 1 , . . . , F ′ n ⊕ F ′′ n α ′′ [b !R := H ]α ′ β ′ [z !R := F ′′ 0 ]β ′′ otherwise. • Let G = F ′ + F ′′ and suppose that bp x,a (F ′ ) = F ′ 0 , F ′ 1 , . . . , F ′ n α ′ β ′ , bp x,a (F ′′ ) = F ′′ 0 , F ′′ 1 , . . . , F ′′ n α ′′ β ′′ , then, bp x,a (F ′ + F ′′ ) := F ′ 0 + F ′′ 0 , F ′ 1 ⊕ F ′′ 1 , . . . , F ′ n ⊕ F ′′ n α ′ α ′′ β ′ β ′′ . Lemma 5.2. Let G be a ground term with fv(G) = x, then bp x,a (G) = O(|G |).
Proposition 5.3 (Soundness of bp). Let G be a ground term whose free variables are given by a sequence x of length n. Then for every r ∈ R n , we have:

bp x,a (G)[a := 1][x := r] (10) n (11)(13)(15) -------------→ O (|G |) ≡ G (r), ∇G(r) .
Proof. Let R = {(10) n (11)(13)(15)}. By induction on G we prove that, supposing bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β, for every vector of real numbers r, we have: We give only the proof of (ii) for G = F ′ [z := F ′′ ], the other cases being similar or trivial. Using the notations from Def. 5.1 we have that the i-th component of the gradient tuple of bp x,a (G) together with the associated substitutions is equal to (we suppose

(i) G 0 β[x := r] R -→ * ≡ G (r); (ii) for all 1 ≤ i ≤ n, for all real number q, G i α β[a := q][x := r] R -→ * ≡ ∂ x i G (r) • q. From (i) and (ii): bp x,a (G)[a := 1][x := r ] ≡ G 0 β[x := r], G 1 α β[a := 1][x := r], . . . , G n α β[a := 1][x := r] R -→ * ≡ G (r), ∇ r G .
F ′ i ⊕ F ′′ i = F ′ i + F ′′ i
, the other cases being simpler):

(F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ [z := F ′′ 0 ]β ′′ [a := q][x := r] ≡ (F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ β ′′ [a := q][x := r][z := F ′′ 0 β ′′ [x := r]] R -→ * (F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)] by (i) ≡ (F ′ i + F ′′ i )α ′′ α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)][b := Hα ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)]] R -→ O (|F ′ |) (F ′ i + F ′′ i )α ′′ α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)][b := ∂ z F ′ (r) • q] by IH ≡ F ′ i α ′ β ′ [a := q][x := r][z := F ′′ (r)] + F ′′ i α ′′ β ′′ [b := ∂ z F ′ (r) • q][x := r] R -→ * ∂ x i F ′ (r) • q + ∂ x i F ′′ (r) • ∂ z F ′ (r) • q by IH (13) ---→ ∂ x i F ′ (r) + ∂ z F ′ (r) • ∂ x i F ′′ (r) • q = ∂ x i G (r) • q
Notice that in order to move to the last line we use the associativity, commutativity and distributivity of + and • over real numbers, not on the corresponding syntactic symbols. In the base cases (variables, functional symbols and numerals), one performs linear substitutions (10) as well as garbage collection (11). In the case of bp x,a (f (x i 1 , . . . ,

x i k ))
one instance of rule ( 15) is needed. □

Notice that the result of the evaluation of bp x,a (G)[a := 1][x := r] is independent from the chosen reduction sequence by Corollary 4.8.

Backpropagation on Higher-Order Programs

Let us now consider the soundness of our transformation ← -D (t) applied to any λ-term of type x !R 1 , . . . , x !R n ⊢ t : R. The proof uses two essential ingredients: first, we remark that the transformation ← -D (t) commutes with any reduction step (Lem. 5.4); second, we prove that ← -D applied to a computational graph G encodes actually all the relevant information of bp x,a (G) (Lem. 5.5). Notice that Lem. 5.5 introduces fresh variables x ′ i (i ≤ n) annotated !R ⊥ , tagging the different components of the gradient of G in a sum. We can then conclude with Th. 5.6 and its Cor. 5.7. Lemma 5.4. Let t be a term of Λ(F ). Then:

(1) if t ≡ t ′ , then ← - D (t) ≡ ← - D (t ′ ), (2) if t ι - → t ′ , for ι any reduction step in Tab. 2, then ← - D (t) X --→ O (1) ← - D (t ′ )
, where X = {ι} for any ι but (13), (15), in these latter cases X = {ι, (8), ( 10), (11)}.

Lemma 5.5. Let G be a ground term with fv(G) = x = x !R 1 , . . . , x !R n and suppose that bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β. Then, there exists J ⊆ {1, . . . , n} such that i J implies G i = 0 and such that

← - D (G)[x ! ← - D (R) := x !R , λa R .x * !R ⊥ a ] - → O (|G |) ≡ G 0 , λa.( j ∈J x * j G j )α β .
Proof. By induction on G. We only consider G = F [z := F ′ ], the other cases being simpler. Let bp x,z,a (F ) = ⟨F 0 , ⟨F 1 , . . . , F n , H ⟩ α⟩ β and bp

x,b (F ′ ) = F ′ 0 , F ′ 1 , . . . , F ′ n α ′ β ′ .
Let us also consider H 0 (the other case being simpler). The term

← - D (G)[x := ⟨x, λa.x * a⟩] is structural equivalent to: ( ← - D (F )[x := x, λa.x * a ])[z := ( ← - D (F ′ )[x := x, λb.x * b ])] - → O (|F ′ |) ≡ ( ← - D (F )[x := x, λa.x * a ])[z := F ′ 0 , λb.( j ∈J ′ x * j F ′ j )α ′ β ′ ] by IH ≡ ( ← - D (F )[x := x, λa.x * a ])[z := F ′ 0 , λb.( j ∈J ′ x * j F ′ j )α ′ ]β ′ (17) ---→ (9) --→ ( ← - D (F )[x := x, λa.x * a ][z := z, z * ])[z := F ′ 0 ][z * := λb.( j ∈J ′ x * j F ′ j )α ′ ]β ′ (16) ---→ ( ← - D (F )[x := x, λa.x * a ][z := z, λa.z * a ])[z := F ′ 0 ][z * := λb.( j ∈J ′ x * j F ′ j )α ′ ]β ′ - → O (|F |) ≡ F 0 , λa.( j ∈J x * j F j + z * H )α β[z := F ′ 0 ][z * := λb.( j ∈J ′ x * j F ′ j )α ′ ]β ′ by IH (10) ---→ (11)
---→ F 0 , λa.( j ∈J

x * j F j + (λb.(

j ∈J ′ x * j F ′ j )α ′ )H )α β[z := F ′ 0 ]β ′ (8)
--→ F 0 , λa.(

j ∈J x * j F j + ( j ∈J ′ x * j F ′ j )α ′ [b := H ])α β[z := F ′ 0 ]β ′ (18)
---→ #J ∩#J ′ F 0 , λa.(

j ∈J ∪J ′ x * j (F j ⊕ F ′ j ))α ′ [b := H ]α β[z := F ′ 0 ]β ′ □
By composing the reductions in Lemma 5.4 and Lemma 5.5, we get:

Theorem 5.6. Let t be a term of

Λ(F ) of type R with fv(t) = x = x !R 1 , . . . , x !R n . For any ground term G such that t ( β -→ ∪ ≡) * G in m β-steps, we have, for a suitable J ⊆ {1, . . . , n}: ← - D (t)[x ! ← - D (R) := x !R , λa R .x * !R ⊥ a ] - → O (m+|G |) ≡ G 0 , λa.( j ∈J x * j G j )α β
where bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β and ∀i J , G i = 0.

Corollary 5.7. Let t be a term of Λ(F ) of type x !R ⊢ t : R, with x = x !R 1 , . . . , x !R n the free variables of t. Let m be the number of β-steps needed to reduce t to a ground term G. For any vector r ∈ R n , let ∇t(r) = д 1 , . . . , д n . Then, for a suitable J ⊆ {1, . . . , n}, with i J , д i = 0, we have:

z * 1[ z, z * := ← - D (t)[x := x, λa.x * a ]][x := r] - → O (m+ |G |) ≡ j ∈J
x * j д j .

One can go further and obtain the tuple of numerals ∇t(r) from j ∈J x * j д j just by: (i) considering ← -D d (t) with d = n (here is the only point where we require a specific choice of d); (ii) replacing each x * i (i ≤ n), morally of type R ⊸ R n , with the corresponding injection λa. ⟨0 . . . , 0, a, 0, . . . , 0⟩; (iii) adding all tuples resulting from the reductions (8),( 10),(11): j ∈J 0, . . . , 0, д j , 0 . . . , 0 -→ * д 1 , . . . , д n . However, this last reduction is quadratic in n, as it performs #J = O(n) sums of vectors of size n without considering that these latter are null but on one coordinate. It would then be preferable to add an ad hoc read-back operation allowing to decode the tuple ∇t(r) out of the tagged sum j ∈J x * j д j more parsimoniously.

Notice that the proof of Cor. 5.7 considers a specific reduction sequence for getting ∇t(r). However, Cor. 4.8 guarantees that the result is independent from the chosen sequence.9 This suggests considering more efficient strategies than that of Cor. 5.7, or even more modular by decomposing the computation along the different components of t. We discuss this last point in the next Section.

AN EXAMPLE: RECURRENT NEURAL NETWORKS

Derivative of a dynamically generated polynomial

Let Nat := (R → R) → R → R be the type of Church natural numbers and consider

G := w • y + x t := λn !Nat .λx .n(λy.G)x .
We have

x !R , y !R , w !R ⊢ G : R and w !R ⊢ t : Nat → R → R. If n encodes n ∈ N,
D (t 2) = (λn.λx .n ← - D (λy.G)x) 2 - → * (λn.λx .n G ′ x) 2 - → λx .2 G ′ x - → λx .G ′ (G ′ x) - → λ x, x * .G ′ (G ′ x, x * ) - → * λ x, x * .G ′ ⟨wx + x, H ⟩ ,
where H := λa.w * (x • a) + x * (w • a) + x * a. The computation continues with

- → * λ x, x * . (w 2 + w + 1)x, λa.w * ((wx + x) • a) + H (w • a) + x * a - → * λ x, x * . (w 2 + w + 1)x, λa.w * ((wx + x) • a) + w * (wx • a) + x * (w 2 • a) + x * (w • a) + x * a - → * λ x, x * . (w 2 + w + 1)x, λa.w * ((2wx + x) • a) + x * ((w 2 + w + 1) • a)
We see that the argument of w * (resp. x * ) is the derivative of f 2 with respect to w (resp. x). This shows how locally free variables are handled correctly. Also observe that G ′ , which is trivial here but could in principle be a very complex function, may be pre-computed independently of n.

Recurrent neural networks

Recurrent neural networks (RNNs) are meant to process inputs that are arbitrary sequences of vectors in R d . They are heavily used for natural language processing: sequences could for instance stand for a word or a sentence, where each vector is a representation of a letter. Such networks iterate over the input to recursively build a desired output. One example is encoding recurrent neural networks, which are used to encode a sequence of vectors x i ∈ R d into an output vector h ∈ R m . For instance, in sentiment analysis [dos Santos and Gatti 2014; [START_REF] Glorot | Domain adaptation for large-scale sentiment classification: A deep learning approach[END_REF][START_REF] Severyn | Twitter sentiment analysis with deep convolutional neural networks[END_REF], it is used to predict if the expressed opinion in a sentence is positive, negative or neutral.

Given a sequence x 1 , ..., x n ∈ R d , the encoding RNN produces intermediate outputs h 1 , ..., h n ∈ R m recursively, by applying a single layer L : R d × R m → R m to both the last output value h i and the next input vector x i+1 :

h i+1 = L(x i+1 , h i ) L(x, h) = σ (E • x + R • h)
where σ is the sigmoid activation function, E is a d × m-matrix called the token embedding matrix, and R is a m × m matrix whose coefficients are called the recurrent weights. h 0 can be initialized to 0. In practice, a RNN ends with a loss function that we are looking to minimize (and which models the problem we want to solve). To keep the example simple, we will not consider it. For the same reason, we suppose that d = m = 1, the general encoding being a direct generalization.

Lists. We represent lists in our language using the Church encoding, which defines a list by its right fold function.

Given A and X be types, the type of lists List(A, X ) is defined as follows:

List(A, X ) := (A → X → X ) → (X → X )
Given a 1 : A, . . . , a n : A, we define the list [a 1 ; ...; a n ] X of type List(A, X ) as:

[a 1 ; ...; a n ] X := λ f !(A→X →X ) .λx !X .f a 1 (f a 2 . . . (f a n x) . . .

)
In what follows we will omit the X annotation if it is clear from context or if it does not matter. Finally, it can be noted that the reverse gradient of the list datatype is given by: ← -

D (List(A, X )) = List( ← - D (A), ← - D (X )) ← - D ([a 1 ; ...; a n ] X ) = [ ← - D (a 1 ), . . . , ← - D (a n )]← - D (X )
Encoding a RNN. Let σ be the function symbol corresponding to the sigmoid function R → R. In dimension 1, the token embedding matrix and the recurrent matrix may be represented as two terms λx !R .(ϵ • x) and λh !R .(ρ • h) respectively, where ϵ !R and ρ !R are two variables. The recurring layer L and the recurring network N is then defined and typed as follows:

L := λx !R .λh !R .σ (ϵ • x + ρ • h) ϵ !R , ρ !R ⊢ L : R → R → R N := λl !List (R,R) .lL0 ϵ !R , ρ !R ⊢ N : List(R, R) → R
The free variables ϵ and ρ are the parameters that we wish to learn via gradient descent.

Backpropagation. For RNNs, the gradient is usually computed using a technique called backpropagation through time [START_REF] Barak | Gradient calculations for dynamic recurrent neural networks: a survey[END_REF][START_REF] Rumelhart | Parallel Distributed Processing[END_REF]]. The method consists in unfolding the RNN (by applying it to an input sequence) and then applying the usual backpropagation over plain vanilla feedforward neural networks. We will now apply our reverse gradient transformation to an example of RNN, and show that backpropagation through time is naturally implemented by the reduction strategy of Theorem 5.6.

Given l = [a 1 ; . . . ; a n ] R , we compute the gradient of Nl with respect to ϵ and ρ. The following proposition exposes the recursive equations expressing the gradient computed thanks to our transformation. They are similar to the equations of backpropagation through time.

Proposition 6.1. We have ∇(Nl)(e, r ) = д n ϵ , д n ρ where д n ϵ and д n ρ are given by the following recurrent equations:

д 0 ϵ = 0 д i+1 ϵ = σ ′ i+1 • (a i+1 + r • д i ϵ ) σ ′ i+1 := ∂ 1 σ (e • a i+1 + r • u i ) д 0 ρ = 0 д i+1 ρ = σ ′ i+1 • (u i+1 + r • д i ρ ) u 0 = 0 u i+1 = σ (e • a i+1 + r • u i )
Proof. By Corollary 5.7, this amounts to computing D = ϵ * д ϵ + ρ * д ρ such that

z * 1[ z, z * := ← - D (
G = F (a 1 , F (a 2 , . . . , F (a n , 0) . . . ))
where L = λx .λh.F with F = ϵ • x + ρ • h. Notice that G is exactly the unfolding of the RNN, which should convince the reader that this strategy implements the backpropagation through time.

We now turn our attention to the effectiveness of the reduction strategy. During the reduction of Nl to G, the λ-abstractions of L must be eliminated. Each of these λ is applied exactly to n different arguments, requiring L to be duplicated n times. In other words, we satisfy the following reduction:

Nl - → * La 1 (La 2 . . . (La n 0) . . . ) - → * G
By point 2 of Lemma 5.4, the exact same reasoning may be applied to ← -D (Nl), in which ← -D (L) must be duplicated n times:

← - D (Nl) - → * ← - D (L) ← - D (a 1 )( ← - D (L) ← - D (a 2 ) . . . ( ← - D (L) ← - D (a n ) ← - D (0)) . . . ) - → * ← - D (G)
A simple observation shows that

← - D (L) is not in β-normal form: ← - D (L) = λx .λh. σ (µ), λa.µ * (∂ 1 σ (µ) • a) [ µ, µ * := ← - D (ϵ • x + ρ • h)] ← - D (ϵ • x + ρ • h)

CONCLUSION AND PERSPECTIVES

On expressiveness. The simply-typed λ-calculus is certainly too restrictive as a programming language, but it does present the main obstacle in defining higher-order backpropagation, namely the native use of higher-order without necessarily going through computational graphs. It also has a minimum of expressiveness for interesting examples to exist, such as inductive types with very basic operations on them (map, fold), as shown in Sect. 6. In fact, our results apply seamlessly to any total language with set-theoretic semantics (e.g. Gödel's System T with arbitrary inductive types), with no need of additional technical ideas. This is already quite broad: no state-of-the-art deep learning architecture we are aware of (convolutional NNs, RNNs, Tree-RNNs, attention-based NNs, transformers. . . ) requires going beyond System T in order to be expressed.

We wish to stress that our current proof does support, unchanged, non-differentiable functions like ReLU. As mentioned at the end of Sect. 2.3, the differentiability hypothesis (⋆) (end of Sect. 4) is essentially cosmetic, in that it allows us to use actual gradients and to avoid the more technical notion of subgradient. In the absence of (⋆), Corollary 5.7 holds for all vectors in the domain of definition of Eq. 25, i.e., wherever the gradient makes sense.

Still in relation with partiality, we argue that our framework can accommodate full recursion (like [START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]), by adding the definition

← - D (Y A ) := Y← - D (A)
to Tab. 3, where Y A : (A → A) → A is the fixpoint operator. As for denotational semantics, one has to consider the cartesian closed category of pointed complete partial order sets and Scott-continuous functions. The type R is interpreted as the flat domain having a bottom element (representing non-termination) and all real numbers as maximal elements. Corollary 5.7 then holds for all vectors on which the program converges and is differentiable. However, this still falls short of being a full programming language like PCF because it lacks conditionals (on R). The issues related to the presence of branchings are mentioned in [START_REF] Plotkin | Some Principles of Differential Programming Languages[END_REF],

where a proposal is suggested for first-order languages.

On complexity. In computational graphs, complexity analysis assumes that computing the value of a single node from its local inputs has unitary cost. In terms of low-level complexity models (e.g. random access machines), this rests on the assumption that searching for the next node to evaluate is a constant-time operation, which is fair because we may suppose the nodes to be linearly ordered compatibly with the dependencies of the graph.

Our analysis shows that, as soon as the programming language is suitably fine-grained, it is consistent to attribute a unitary cost to a single evaluation step (i.e., a rewriting rule from Tab. 2). However, in general programming languages, it is unclear whether the search for the next redex has constant cost. Albeit recent work [Accattoli and Barras 2017, Corollary 9.2] shows that this is the case in call-by-name evaluation of closed programs, a detailed analysis for our language is currently missing and we defer it to future work.

On higher-type derivatives. A major endeavor relating functional primitives with differential operators is differential linear logic [START_REF] Ehrhard | An introduction to differential linear logic: proof-nets, models and antiderivatives[END_REF]] and its associated differential λ-calculus [START_REF] Ehrhard | The differential lambda-calculus[END_REF]. Indeed, the initial motivation of our work was to express the backpropagation algorithm in the differential λ-calculus. Although this is possible, we realized that it required discarding the most important programming primitive of the differential λ-calculus, namely the derivative operator D on higher-order types. This is probably a good place to point out a crucial feature of our program transformation ← -D . Take a term λx .t(ux) of type R → R and such that ⊢ t : A → R and ⊢ u : R → A. We have two different ways of computing the derivative of t(ux) with respect to x R : either by our transformation ← -D (λx .t(ux)) or by Ehrhard's derivative operator D(λx .t(ux)). Both ways are purely functional, but our operator decomposes as natural and intriguing question to understand whether the ability to compute derivatives at higher types (rather than just over R) can play a role in developing new machine learning models.

← - D (λx .t(ux)) = β λ ⟨x, a⟩ . ← - D (t)( ← - D (
The explicit substitution rule is similar to the previous one, using (10) or (11). If t = f (t 1 , . . . , t n ), then we have by IH that each t i is a closed value, so a numeral, so we can apply the numerical rule associated with f . If t is a numeral, then it is a value. If t = t 1 + t 2 , then by IH each t i is a value so we may apply either ( 14) or (13). □ Proposition 4.5 For every term t, and every set X ⊆ βℓ, there exists a X -normal form u such that t X --→ * u.

Proof (Sketch). In case (10) X , it is trivial to find a notion of weight w(t) ∈ N such that t X --→ t ′ implies w(t) > w(t ′ ). Otherwise, let k be the greatest size of the type of a subterm of t, and let W (t) be the sequence (n 1 , . . . , n k , n k +1 ) where n k+1 = w(t) and for every 1 ≤ i ≤ k, n i is the number of redexes of type ( 8), ( 9), ( 10) in t such that the term in the explicit substitution has a type of size k + 1 -i (so that, as i grows bigger, n i accounts for smaller and smaller types). It is not hard to show that, if t is not X -normal, then t X --→ t ′ for some t ′ such that W (t) > W (t ′ ) in the lexicographical order. □ Lemma 4.6 For any t (10) n (11)(13)( 15)

-------------→ * u, the number of reduction steps in the sequence is O(|t |).

Proof. Let ∥t ∥ be the number of symbols in t where each variable is counted twice. The claim then follows by remarking that any rule ( 11), ( 13), ( 15) makes ∥t ∥ strictly decrease, and the same for (10) whenever the substituted term is a numeral. □ Table 4 recalls the denotational model induced by the cartesian closed structure of the category of sets and functions.

The proof of the soundness property is completely standard.

A.2 Proofs of Section 5

Lemma 5.2. Let G be a ground term whose free variables are given by a sequence x of length n. Then, By induction on the definition of bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β, one proves that:

(i) |G 0 β | = O(|G |); (ii) for all 0 < i ≤ n, |G i | = O(k(# x i G + 1)), with # x i G is the number of occurrences of x i in G; (iii) for all 0 < i ≤ n, |G i α | = O(|G |). Notice that the items (i)-(iii) give that bp x,a (G) = n i=1 |G i | + |G 0 β | + |α | = O(# x 1 G + • • • + # x n G) + 2 |G |) = O(|G |)). We detail only the case G = F [z := F ′ ],
taking the notation of the definition of bp( ). Item (i) follows trivially from the induction hypothesis. Item (ii) is a consequence of the induction hypothesis and the fact that #

x i (G) = # x i (F ) + # x i (F ′ ).
Concerning item (iii), we have that:

(F i ⊕ F ′ i )α ′ [b := H ]α = |F i α | + F ′ i α ′ + 1 + |H | = O(|F | + |F ′ |) = O(|G |). The factor k in (ii) is due to the case G = f (x i 1 , . . . , x i ′ k ). □ R := R A × B := A × B A → B := set of functions from A to B R ⊥ d := set of linear maps from R to R d x !A 1
Proposition 5.3. Let G be a ground term whose free variables are given by a sequence x of length n. Then for every r ∈ R n , we have: bp

x,a (G)[a := 1][x := r] (10) n (11)(13)(15) -------------→ O (|G |) ≡ G (r), ∇G(r) .
Proof. By induction on G we prove that, if bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β, then, for every vector of real numbers r, we have:

(i) G 0 β[x := r] (10) n (11)(13)(15) -------------→ * ≡ G (r); (ii) for all 1 ≤ i ≤ n, for all real number q, G i α β[a := q][x := r] (10) n (11)(13)(15) -------------→ * ≡ ∂ x i G (r) • q;
From (i) and (ii) then follows: -------------→ * ∂

x i F ′ (r) • q + ∂ x i F ′′ (r) • ∂ z F ′ (r) • q by IH (13) ---→ ∂ x i F ′ (r) • q + ∂ x i F ′′ (r) • ∂ z F ′ (r) • q = ∂ x i F ′ (r) + ∂ z F ′ (r) • ∂ x i F ′′ (r) • q = ∂ x i G (r) • q
Notice that in order to move to the last line we use the associativity, commutativity and distributivity of + and • over real numbers, not on the corresponding syntactic symbols.

Let us consider also the case G = ⟨F ′ , F ′′ ⟩. Using the notations from the definition:

G 0 β[x := r] = F ′ 0 , F ′′ 0 β ′ β ′′ [x := r] ≡ F ′ 0 β ′ [x := r], F ′′ 0 β ′′ [x := r]
(10) n (11)(13)(15)

-------------→ * ≡ F ′ (r), F ′′ (r) by IH

= G (r)

The case (ii) is similar (just a sum instead of a pair).

In the base cases (variables, functional symbols and numerals), one performs linear substitutions (10) as well as garbage collection (11). In the case of bp x,a (f (x i 1 , . . . , x i k )) one instance of the rule ( 15) is needed. □ Lemma A.1. Let t be a term of Λ(F ). We have:

← - D (t) = O(|t |).
Proof. By just inspecting the definition in Table 3. □ Lemma 5.4. Let t be a term of Λ(F ). Then:

(1) if t ≡ t ′ , then

← - D (t) ≡ ← - D (t ′ ),
(2) if t ι -→ t ′ , for ι any reduction step in Tab. 2, then ← -D (t) X --→ O (1) ← -D (t ′ ), where X = {ι} for any ι but (13), ( 15), in these latter cases X = {ι, (8), ( 10), (11)}.

Proof. Item 1 is immediate by case inspection. Concerning item 2, one proves by induction on t that: Then, suppose t ι -→ t ′ , notice that ι (18) since t has no R ⊥ -variable. If ι = (10), then the statement follows from (•). All other cases are immediate. We detail just the case of an instance of a numerical rule (13), ( 15), e.g. r + q -→ r + q. We have that ← -D (r + q) is:

x + y, λa.(x ′ a + y ′ a) [ x, x ′ := r , λa.0 ][ y, y ′ := q, λa.0 ] (8)(10)(11)

--------→ O (1) r + q, λa.(0 + 0)

---→ 2 r + q, λa.0 □ Lemma 5.5. Let G be a ground term of type x !R ⊢ G : R with x = x !R 1 , . . . , x !R n the free variables of G and let bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β. For a suitable J ⊆ {1, . . . , n}, such that if i J , then G i = 0, we have:

← - D (G)[x ! ← - D (R) := x !R , λa R .x ′!R ⊥ a ] - → O (|G |) ≡ G 0 , λa.( j ∈J
x ′ j G j )α β .

Proof. By induction on G. We will freely use structural equivalence where necessary, knowing that Proposition x ′ j F ′ j )α ′ ]β ′

Now we can apply the induction hypothesis on ← -D (F ) and we split in two sub-cases, depending whether the variable z ′ appears in the tuple associated with ← -D (F ) or not. In the first case, we have that the above term reduces by IH to:

- → O (|F |) ≡ F 0 , λa.( j ∈J x ′ j F j + z ′ H )α β[z := F ′ 0 ][z ′ := λb.( j ∈J ′ x ′ j F ′ j )α ′ ]β ′ by IH (10)
---→

---→ F 0 , λa.( j ∈J

x ′ j F j + (λb.(

j ∈J ′ x ′ j F ′ j )α ′ )H )α β[z := F ′ 0 ]β ′ (8)
--→ F 0 , λa.( j ∈J

x ′ j F j + (

j ∈J ′ x ′ j F ′ j )α ′ [b := H ])α β[z := F ′ 0 ]β ′ (18)
---→ #J ∩#J ′ F 0 , λa.(

j ∈J ∪J ′ x ′ j (F j ⊕ F ′ j ))α ′ [b := H ]α β[z := F ′ 0 ]β ′
Notice that the above instance of (10) replaces exactly one occurrence of z ′ . The last term satisfies the statement of the lemma.

In the case the variable z ′ does not appear in the tuple associated with ← -D (F ), this means that H = 0. We then have:

β ηℓ ---→ O (|F |) ≡ F 0 , λa.( j ∈J x ′ j F j )α β[z := F ′ 0 ][z ′ := λb.( j ∈J ′ x ′ j F ′ j )α ′ ]β ′ by IH - → F 0 , λa.( j ∈J x ′ j F j )α β[z := F ′ 0 ]β ′ by (11)
The last term satisfies the statement of the lemma.

Let G = f (x i 1 , . . . , x i k ) for a suitable subset x i 1 , . . . , x i k ⊆ x. In this case we have:

← -D (G)[x := x, λa.x ′ a ] = f (y) , λa. x ′ j F j )α]β[y ′ 2 := λa.(

j ′ ∈J ′ x ′ j ′ F ′ j ′ )α ′ ]β ′ by IH (10)(8) (11)
--------→ O (1) F 0 + F ′ 0 , λa.((

j ∈J x ′ j F j )α + ( j ′ ∈J ′ x ′ j ′ F ′ j ′ )α ′ ) ββ ′ ≡ F 0 + F ′ 0 , λa.(( j ∈J
x ′ j F j ) + (

j ′ ∈J ′ x ′ j ′ F ′ j ′ ))αα ′ ββ ′ (18)
---→ #J ∩#J ′ F 0 + F ′ 0 , λa.(

j ∈J ∪J ′ x ′ j F j ⊕ F ′ j ′ )αα ′ ββ ′
Notice that all the above instances of (10) replaces exactly one occurrence of a variable. The last term satisfies the statement of the lemma. □ Theorem 5.6. Let t be a term of Λ(F ) of type x !R ⊢ t : R, with x = x !R 1 , . . . , x !R n the free variables of t. For any ground term G such that t (

  with n inputs and m outputs invoking |G | operations (i.e., nodes), forward mode computes the Jacobian of G in O(n |G |) operations, while reverse mode requires O(m |G |) operations. In deep learning, as the number of layers increases, n

  by-value is the same as executing the forward mode AD algorithm on G. Moreover, the definition of -→ D is fully compositional. 3 Indeed, we merely reformulated the transformation --→ D x introduced in[START_REF] Wang | Demystifying differentiable programming: shift/reset the penultimate backpropagator[END_REF]: the computation of the proposition corresponds to that of ---→ D x j (G).

(

  18), allowing us to evaluate ← -D (G) with the same efficiency as bp x 1 ,x 2 ,a (G). The linear factoring f t + f u -→ f (t + u) would be semantically unsound in general if f : ¬R = R → R, which is why we must explicitly track the linearity of negations. So we have a compositional transformation ← -D which takes a computational graph G with n inputs x 1 , . . . , x n and returns a program x 1 : R × R ⊥ , . . . , x n : R × R ⊥ ⊢ ← -D (G) : R × R ⊥ in the simply-typed λ-calculus augmented with linear negation, such that ← -D (G) evaluates to (essentially) bp x 1 , ...,x n ,a (G). Actually, thanks to another nice observation of [Wang et al. 2019], we can do much more: we can extend ← -D for free to the whole simply-typed λ-calculus, just letting ← -D (fun x ->t ) := fun x -> and, whenever x 1 : R, . . . , x n : R ⊢ t : R, we have that ← -D (t)

X

  --→ u. If the set X is a singleton {ι}, we simply writeι -→. If X = β ∪ η ∪ ℓ, we write just -→. If k ∈ N, X --→ k denotes a sequence of length k ofX denotes the context closure of all reduction relations except (18).

  Let us underline that, in the last line of Tab. 3b, the index d of ← -D d is totally independent from the arity k of the function symbol f and from the indexes i ≤ k of the variables x * i tagging the sum introduced by ← -D d . The fact that d

  Prop. 4.3 allows to postpone all structural equivalences at the end. From Lem. 4.6 and bp x,a (G) = O(|G |) (Lem. 5.2) the length of the reduction is O(|G |).

  the term t n dynamically generates the function λx .f n (w, x) with f n (w, x) := (w n + w n-1 + • • • + w + 1) • x (we take some liberty in simplifying arithmetic expressions for the sake of readability). Notice that the free variable x of λy.G, a term to be iterated at will, is captured later in t. We have ← -D (λy.G) -→ * λ y, y * . G, λa.w * (y • a) + y * (w • a) + x * a =: G ′ , where we used the shorthand λ⟨y, y * ⟩.u := λp.u[⟨y, y * ⟩ := p]. Hence, when e.g. n = 2, ← -

  Nl)[ϵ := ϵ, λa.ϵ * a ][ρ := ρ, λa.ρ * a ]][ϵ := e][ρ := r ] -→ * D This is done by induction over the size of the list l. □ The strategy from Theorem 5.6 first computes ← -D (Nl), then reduces it to ← -D (G) such that Nl -→ * G, with G a ground term, i.e., only containing subterms of ground type. Here:

  being a pair, we have ← -D (L) β -→ k u for some β-normal M. Moreover, it can be shown that when the dimensions d and m are > 1, k = O((d + m) × m). In our case, we see that each component of the sum in← -D (ϵ • x + ρ • h)will be responsible for at least one substitution. Finally, since ← -D (L) is duplicated n times in the original backpropagation through time strategy, by using the strategy where ← -D (L) is reduced to M before being substituted, we obtain a gain of O((d +m)mn) reduction steps. This is significant as the number of learning parameters (here (d +m)m) can be very large in typical neural networks. For instance, a recent model called GPT-2[START_REF] Radford | Language models are unsupervised multitask learners[END_REF]] has 1.5 billion parameters.

  bp x,a (G) = O(|G |). Proof. We actually prove a slight stronger claim: let G be a ground term whose free variables are in a sequence x of length n. Then, bp x,a (G) = O(n + |G |). Notice that then whenever x contains exactly the free variables of G, n ≤ |G | and so O(n + |G |) = O(|G |).

  bp x,a (G)[a := 1][x := r ] ≡ G 0 β[x := r], G 1 α β[a := 1][x := r], . . . , G n α β[a := 1][x := r] (10) n (11)(13)(15)-------------→ * ≡ G (r), ∇G(r) .

  Proposition 4.3 allows to postpone all structural equivalences at the end. From Lemma 4.6 and the fact that bp x,a (G) = O(|G |) (Lemma 5.2) the length of the reduction is O(|G |)The non-trivial case for (i) and(ii) is G = F ′ [z := F ′′ ].Using the notations from the definition:G 0 β[x := r] = F ′ 0 β ′ [z := F ′′ 0 ]β ′′ [x := r] ≡ F ′ 0 β ′ [x := r][z := F ′′ 0 β ′′ [x := r]] (10) n (11)(13)(15) -------------→ * F ′ 0 β ′ [x := r][z := F ′′ (r)] by IH (10) n (11)(13)(15) -------------→ * F ′ (r, F ′′ (r)) = G (r) by IHAs for (ii), we have that the i-th component of the gradient tuple of bp x,a (G) together with the associated substitutions is equal to (we supposeF ′ i ⊕ F ′′ i = F ′ i + F ′′ i, the other cases being simpler):(F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ [z := F ′′ 0 ]β ′′ [a := q][x := r] ≡ (F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ β ′′ [a := q][x := r][z := F ′′ 0 β ′′ [x := r]] (10) n (11)(13)(15) -------------→ * (F ′ i + F ′′ i )α ′′ [b := H ]α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)] by (i) ≡ (F ′ i + F ′′ i )α ′′ α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)][b := Hα ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)]] (10) n (11)(13)(15) -------------→ O (|F |) (F ′ i + F ′′ i )α ′′ α ′ β ′ β ′′ [a := q][x := r][z := F ′′ (r)][b := ∂ z F ′ (r) • q] by IH ≡ F ′ i α ′ β ′ [a := q][x := r][z := F ′′ (r)] + F ′′ i α ′′ β ′′ [b := ∂ z F ′ (r) • q][x := r](10) n (11)(13)(15)

(

  4.3allows us to postpone all structural equivalences to the end without affecting the number of reduction steps. The non-trivial cases are detailed below. Remark that the length of each reduction is O(← -D (G) ) which is equal to O(|G |) by Lemma A.1. Let G = F [z := F ′ ] and suppose that bp x,z,a (F ) = ⟨F 0 , ⟨F 1 , . . . , F n , H ⟩ α⟩ β, bp x,b (F ′ ) = F ′ 0 , F ′ 1 , . . . , F ′ n α ′ β ′ .We have that← -D (G)[x := ⟨x, λa.x ′ a⟩] is structural equivalent to : ( ← -D (F )[x := x, λa.x ′ a ])[z := ( ← -D (F ′ )[x := x, λb.x ′ b ])] -→ O (|F ′ |) ≡ ( ← -D (F )[x := x, λa.x ′ a ])[z := F ′ 0 , λb.( [x := x, λa.x ′ a ])[z := F ′ 0 , λb.( [x := x, λa.x ′ a ][z := z, z ′ ])[z := F ′ 0 ][z ′ := λb.()[x := x, λa.x ′ a ][z := z, λa.z ′ a ])[z := F ′ 0 ][z ′ := λb.( j ∈J ′

  f (y) • a [ y, y ′ := x][x := x, λa.x ′ a ] (9)(10) -----→ O (n) f (y) , λa. k j=1 y ′ j ∂ j f (y) • a [y := x][y ′ := λa.x ′ a] ≡ f (y) , λa. k j=1 y ′ j ∂ j f (y) • a [y := x][y := x][y ′ := λa.x ′ a] (8)(10) -----→ O (k ) (11) ---→ O (n) f (x i 1 , . . . , x i k ) , λa. k j=1 x ′ i j ∂ j f (x i 1 , . . . , x i k ) • a Let G = F + F ′ and suppose that bp x,a (F ) = ⟨F 0 , ⟨F 1 , . . . , F n ⟩ α⟩ β, bp x,a (F ′ ) = F ′ 0 , F ′ 1 , . . . , F ′ n α ′ β ′ .We have that← -D (G)[x := ⟨x, λa.x ′ a⟩] is equal to:y 1 + y 2 , λa.(y ′ 1 a + y ′ 2 a) [ y 1 , y ′ 1 := ← -D (F )][ y 2 , y ′ 2 := ← -D (F ′ )][x := x, λa.x ′ a ] β ηℓ ---→ O (|G |) ≡ F 0 + F ′ 0 , λa.(y ′ 1 a + y ′ 2 a) [y ′ 1 := λa.( j ∈J

β

  -→ ∪ ≡) * G in m β-steps, we have, for a suitable J ⊆ {1, . . . , n}:← -D (t)[x := x, λa.x ′ a ] -→ O (m+|G |) ≡ G 0 , λa.( j ∈J x ′ j G j )α β where bp x,a (G) = ⟨G 0 , ⟨G 1 , . . . , G n ⟩ α⟩ β and ∀i J , G i = 0.Proof. Let us suppose t ( β -→ ∪ ≡) * G in m steps. By Lemma 5.4, we have that:← -D (t)[x := x, λa.x ′ a ] ( β -→ ∪ ≡) * ← -D (G)[x := x, λa.x ′ a ]

  ) is sound. The interpretation of a judgment Γ ⊢ t : A is a function t Γ from the cartesian product Γ of the denotations of the types in Γ to A . The interpretation of a judgment Γ ⊢ z t : R d , is given as a function t Γ;z associating with every g ∈ Γ a linear map from R to R d . The definition is by induction on t and completely standard (explicit substitution is functional composition). We omit the superscript Γ (or Γ; z) if irrelevant. This interpretation supposes to have associated each function symbol in F with a suitable map over real numbers. By a standard reasoning, one can prove that: Proposition 4.7 (Semantic soundness). Let Γ ⊢ (z) t : A, then t -→ u or t ≡ u, gives t = u .The semantic soundness gives as a by-product a light form of confluence on Euclidean types 8 .Corollary 4.8. Let ⊢ t : R d and v and v ′ be β-normal forms s.t. t -→

* v and t -→

Table 3 .

 3 The reverse gradient ← -D

d relative to R d , for an arbitrary natural number d .

  u) ⟨x, a⟩), while Ehrhard's follows the chain rule (see Sect. 2.3), giving D(λx .t(ux)) = β λ ⟨x, a⟩ .D(t) ⟨ux, D(u) ⟨x, a⟩⟩. 10 It is clear that this latter expression is inefficient with respect to backpropagation because of the duplication of the term u (see also Sect. 2). D and the operator D of the differential λ-calculus are extensionally different. This is immediately seen on "pure" λ-terms (with no function symbols): when A is higher order, there are pure terms t : A → R with non-trivial derivative, which is always computed by D(t), whereas ← -D (t) = t (modulo a type change), which cannot compute the derivative. 11 Indeed, observe that our soundness statements (Theorem 5.6, Corollary 5.7) only hold for terms of ground type. At present, we seem to have no use for the extra generality provided by the differential λ-calculus, but it is a

	← -
	So our

If a category C has biproducts, then 1 0 (the terminal object is also initial). If C is cartesian closed, then its product × is a left adjoint and therefore preserves colimits (initial objects in particular), so A A × 1 A × 0 0 for every object A of C.

i ∈I t i , 14 Γ ⊢ z z : R Γ, x !A ⊢ x : A Γ ⊢ (z) t : A Γ ⊢ (z) u : B Γ ⊢ (z) ⟨t, u ⟩ : A × B Γ ⊢ u : A × B Γ, x !A , y !B ⊢ (z) t : C Γ ⊢ (z) t [ x !A , y !B := u] : C Γ, x !A ⊢ t : B Γ ⊢ λx !A .t : A → B Γ ⊢ t : A → B Γ ⊢ u : A Γ ⊢ tu : B Γ ⊢ z t : R d Γ ⊢ λz R .t : R ⊥ d Γ ⊢ t : R ⊥ d Γ ⊢ (z) u : R Γ ⊢ (z) tu : R d Γ ⊢ u : A Γ, x !A ⊢ (z) t : C Γ ⊢ (z) t [x !A := u] : C Γ ⊢ (z ′ ) u : R Γ ⊢ z t : R d Γ ⊢ (z ′ ) t [z R := u] : R d Γ ⊢ t 1 : R . . . Γ ⊢ t k : R Γ ⊢ f (t 1 , . . . , t k ) : R r ∈ R Γ ⊢ r : R Γ ⊢ (z) t : R Γ ⊢ u : R Γ ⊢ (z) t • u : R Γ ⊢ t : R Γ ⊢ (z) u : R

In this paper we focus on exactly what is required to express the backpropagation algorithm in the λ-calculus, avoiding a full linear logic typing assignment and just tracking the linearity of a single variable of type R in judgments typing a term with a Euclidean type R d .

More sophisticated notions of confluence modulo an extension of ≡ hold for the reductions in Tab. 2, but we avoid to discuss this point here because inessential for our purposes.

This is in fact the case if one admits the injections λa . ⟨0 . . . , 0, a, 0, . . . , 0⟩ of type R ⊥ , as discussed in the previous note.

Here we use some syntactic sugar in order to avoid notational bureaucracy. For example, λ ⟨x, a ⟩ . . . stands for λy . . . . [⟨x, a ⟩ := y] in our language. Also, the differential λ-calculus of[START_REF] Ehrhard | The differential lambda-calculus[END_REF] does not have pairs and so λ ⟨x, a ⟩ is curried into λx .λa.

This does not contradict Corollary 5.7 because, when A is of order zero (i.e., A = R n ), any pure t must be a projection and the above equality is correct.

, . . . , x !A n n := A 1 × • • • × A n (a) interpretation of types Γ ⊢ t : A : function from Γ to A Γ ⊢ z t : R d : function from Γ to R ⊥ d Γ ⊢ x : A (g) := д x (27) Γ ⊢ z z : R (g) := r → r (28) Γ ⊢ f (t) : R (g) := f ( Γ ⊢ t : R (g)) (29) Γ ⊢ r : R (g) := r (30) Γ ⊢ λx !A .t : A → B (g) := d → Γ, x !A ⊢ t : B (g, d)(31)Γ ⊢ tu : B Γ (g) := Γ ⊢ t : A → B (g)( Γ ⊢ u : A (g)) (32) Γ ⊢ λz R .t : R ⊥ d (g) := Γ ⊢ z t : R d (g) (33) Γ ⊢ z tu : R d (g) := a → Γ ⊢ t : R d (g)( Γ ⊢ z u : R (g)(a))(34)Γ ⊢ z ⟨t, u⟩ : R d × R d ′ Γ (g) := a → ( Γ ⊢ z t : R d (g)(a), Γ ⊢ z u : R d ′ (g)(a))(35)Γ ⊢ z t[x !A := u] : R d (g) := a → Γ, x !A ⊢ t : R d (g, Γ ⊢ u :: R (g))(a) (36)Γ ⊢ z 1 t[z R

:= u] : R d (g) := a → Γ ⊢ z 2 t : R d (g)( Γ ⊢ z 2 u : R (g)(a)) (37) Γ ⊢ z t[ x !A , y !B := u] : R d (g) := a → Γ, x !A , y !B ⊢ z t : R d (g, Γ ⊢ u : A × B (g))(a)(38)Γ ⊢ z t • u : R (g) := a → Γ ⊢ z t : R (g)(a) • Γ ⊢ u : R (g)

ACKNOWLEDGMENTS

We would like to thank T. Ehrhard, C. Fouqueré, M. Gaboardi, B.A. Pearlmutter, Y. Regis-Gianas, J.M. Siskind, C. Tasson and the anonymous reviewers for useful comments and discussions.

A APPENDIX

A.1 Proofs of Section 4 Proposition 4.1 If t -→ u or t ≡ u and Γ ⊢ t : A, then Γ ⊢ u : A.

Proof (Sketch). Let C{t } -→ C{u}(resp. C{t } ≡ C{u}), where t -→ u (resp. t ≡ u) is a reduction step (resp. a structural equivalence) in Tab. 2. The proof is a standard induction on the evaluation context C. The base case is by inspection of all cases in Tab. 2, using the following properties:

• (substitution) if Γ, x A ⊢ t : B and Γ ⊢ u : A then Γ ⊢ t {u/x } : B;

• (splitting) if Γ, x A ⊢ t : B then for any term t {y/x } obtained from t by renaming some (possibly none) occurrences of x to a fresh variable y, we have Γ, x A , y A ⊢ t {y/x } : B.

The two properties are achieved by induction on t. Notice that the splitting property implies the so-called "weakening" lemma: Γ ⊢ t : B implies Γ, x A ⊢ t : B for a fresh variable x. □ Lemma 4.2 Let ι be any reduction rule and let

Proof (Sketch). The proof goes along the same lines as [START_REF] Accattoli | Distilling Abstract Machines[END_REF]]. We denote by ≡ 1 the symmetric and context closure of the relation induced by the rules ( 19)-( 23), i.e., ≡ is the reflexive-transitive closure of ≡ 1 . We then prove the following lemmata.

(1) If v ≡ 1 t and v is a value, then t is also a value.

(Remark that no structural equivalence rule allows for commuting an explicit substitution with a λ-abstraction.

In the case v = ⟨v 1 , v 2 ⟩, then the only possibility is u = ⟨u 1 , u 2 ⟩ α (rules ( 22), ( 23)), but in this case we would have v i = u i α, for i = 1 or i = 2, contradicting the hypothesis of v i to be a value).

(2) If C{x } ≡ 1 u, then u = C ′ {xα }, for some context C ′ and explicit substitution α.

(By induction on C and case inspection).

(3) If t ≡ 1 t ′ and t ι -→ u, for some reduction step applied to the root of t, then there exists u ′ ≡ u such that t ′ ι -→ u ′ .

(By case inspection. Notice that in case t

---→ t and t[x

(1) to assure that u is a value, in particular no explicit substitution has been lifted from inside v).

(

(By 2 and 3).

The general statement then follows by induction on the length of a sequence t

for the induction step. □ Proposition 4.3 Let X be any subset of the reduction rules in Tab. 2 (including the variant (10) n ) and let t (

Proof. By induction on k, using Lemma 4.2. □ Proposition 4.4. Given a closed term t of type A, if t is a β-normal form, then it is a value.

Proof. By induction on a derivation ⊢ t : A. If the last rule is an → introduction, then t is a value. If it is an → elimination, then t = u 1 u 2 and ⊢ u 1 : B → A or ⊢ u 1 : R ⊥ . By IH u 1 is then a value, hence of the form λx .u ′ and therefore u 1 u 2 is a redex (8). If the last rule is a × introduction, we apply trivially the IH. If it is a × elimination, then t = t ′ [⟨x, y⟩ := u] and ⊢ u : A × B. By IH u is a value, and hence it is of the form u = ⟨v 1 , v 2 ⟩, therefore t is a (9) redex. (39) Γ ⊢ z 0 : R (g) := a → 0 (40)

(b) interpretation of terms. In the cases admitting two version (with or without linear variables), we considered only the case with linear variables, the other being simpler. The reader may check that the linearity (i.e. commutation with sums and scalar products) of the parameter a associated with a linear variable is respected.

Table 4. The denotational model induced by the category of sets and functions.

Notice that the above lemma would fail if we replace the meta-operator ⊕ with the sum constructor + in the definition of bp x,a (G). In that case we would have the size of bp x,a (x

copies the n -1 sums for every component of the gradient tuple.

in O(m) steps. By Lemma 5.5, we have that:

x ′ j G j )α β this latter term satisfying the statement of the theorem. We compose the reductions and apply Proposition 4.3 to conclude. □ Corollary 5.7. Let t be a term of Λ(F ) of type x !R ⊢ t : R, with x = x !R 1 , . . . , x !R n the free variables of t. Let m be the number of β-steps needed to reduce t to a ground term G. For any vector r ∈ R n , let ∇t(r) = д 1 , . . . , д n . Then, for a suitable J ⊆ {1, . . . , n}, with i J , д i = 0, we have:

Proof. Apply Theorem 5.6 to ← -D (t)[x := ⟨x, λa.x ′ a⟩] and then Proposition 5.3.

where

Then the result follows by applying recursively Lemma A.2. □ Proposition 6.1. We have ∇(Nl)(e, r ) = д n ϵ , д n ρ where д n ϵ and д n ρ are given by the following recurrent equations: We get the result by direct application of Lemma A.3. □