Patent Publication Number: US-2022237258-A1

Title: System and Method for Automatic Differentiation of Higher-Order Functions

Description:
This application claims the benefit of priority of U.S. provisional application Ser. No. 62/840,387, filed on Apr. 29, 2019 the disclosure of which is herein incorporated by reference in its entirety. 
    
    
     GOVERNMENT LICENSE RIGHTS 
     This invention was made with government support under contract numbers IIS51522954 and U.S. Pat. No. 1,734,938 awarded by the National Science Foundation. The government has certain rights in the invention. 
    
    
     FIELD 
     The device and method disclosed in this document relates to computer processing systems and, more particularly, to a system and method for automatic differentiation of higher-order functions. 
     BACKGROUND 
     Unless otherwise indicated herein, the materials described in this section are not prior art to the claims in this application and are not admitted to the prior art by inclusion in this section. 
     The classical univariate derivative of a function ƒ: →  is a function ƒ′: → . Multivariate or vector calculus extends the notion of derivative to functions whose domains and/or ranges are aggregates, that is vectors, introducing notions like gradients, Jacobians, and Hessians. Differential geometry further extends the notion of derivatives to functions whose domains and/or ranges are—or can contain—functions. 
     Automatic differentiation (AD) is a collection of methods for computing the derivative of a function at a point when the function is expressed as a computer program. These techniques, once pursued mainly by a small quiet academic community, have recently moved to the forefront of deep learning, where more expressive languages can spawn new industries, efficiency improvements can save billions of dollars, and errors can have far-reaching consequences. 
     From its earliest days, AD has supported functions whose domains and/or ranges are aggregates. However, there is currently interest from application programmers (machine learning in particular) in applying AD to higher-order functions. Classical AD systems, such as ADIFOR, TAPENADE, were implemented for first-order languages like FORTRAN, C, and C++. However, such classical AD systems do not expose a derivative-taking operator as a higher-order function, let alone one that can take derivatives of higher-order functions. More recent AD systems, such as MYIA, and TANGENT, as well as the HASKELL AD package available on Cabal, the “Beautiful Differentiation” system, and the “Compiling to Categories” system, have been implemented for higher-order languages like SCHEME, ML, HASKELL, F , PYTHON, LUA, and JULIA. All these recent systems expose a derivative-taking operator as a higher-order function. However, they do not support taking derivatives of higher-order functions. Accordingly, what is needed is an AD system that exposes a derivative-taking operator as a higher-order function and supports taking derivatives of higher-order functions, and does so in a manner that reliably yields correct results. 
     SUMMARY 
     A method for automatic differentiation of a function defined by program code is disclosed. The method comprises storing, in at least one memory, an automatic differentiation program having a derivative operator program construct that implements a mathematical derivative operator. The method comprises receiving, with at least one processor, first program code defining a first mathematical function, the first mathematical function being a higher-order function which takes a first argument as input. The method comprises determining, by executing the automatic differentiation program with the at least one processor, a second mathematical function, the second mathematical function being a derivative of the first mathematical function, the determining of the second mathematical function including multiple executions of the derivative operator program construct, each execution of the derivative operator program construct being distinguished from one another using unique distinguishing features. 
     A system for automatic differentiation of a function defined by program code is disclosed. The system comprises at least one memory configured to store program instructions including an automatic differentiation program having a derivative operator program construct that implements a mathematical derivative operator. The system comprises at least one processor configured to execute the program instructions stored on the at least one memory to: receive first program code defining a first mathematical function, the first mathematical function being a higher-order function which takes a first argument as input; and determine a second mathematical function, the second mathematical function being a derivative of the first mathematical function, the determination of the second mathematical function including multiple executions of the derivative operator program construct, each execution of the derivative operator program construct being distinguished from one another using unique distinguishing features. 
     A non-transitory computer readable medium for automatic differentiation of a function defined by program code is disclosed. The computer-readable medium stores program instructions including an automatic differentiation program having a derivative operator program construct that implements a mathematical derivative operator, the program instructions, when executed by at least one processor, cause the at least one processor to: receive first program code defining a first mathematical function, the first mathematical function being a higher-order function which takes a first argument as input; and determine a second mathematical function, the second mathematical function being a derivative of the first mathematical function, the determination of the second mathematical function including multiple executions of the derivative operator program construct, each execution of the derivative operator program construct being distinguished from one another using unique distinguishing features. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The foregoing aspects and other features of the system and method are explained in the following description, taken in connection with the accompanying drawings. 
         FIG. 1  shows a block diagram of an exemplary embodiment of an automatic differentiation system. 
         FIG. 2  shows a flow diagram for a method for operating the automatic differentiation system. 
         FIGS. 3A-3B  show program code for a first exemplary implementation of the automatic differentiation program. 
         FIGS. 4A-4C  show program code for a second exemplary implementation of the automatic differentiation program. 
     
    
    
     DETAILED DESCRIPTION 
     For the purposes of promoting an understanding of the principles of the disclosure, reference will now be made to the embodiments illustrated in the drawings and described in the following written specification. It is understood that no limitation to the scope of the disclosure is thereby intended. It is further understood that the present disclosure includes any alterations and modifications to the illustrated embodiments and includes further applications of the principles of the disclosure as would normally occur to one skilled in the art which this disclosure pertains. 
     Introduction to Automatic Differentiation of Higher-Order Functions and its Challenges 
     As discussed in greater detail below, the disclosure provides two methods that that enable automatic differentiation of higher-order functions whose ranges and/or domains are functions. Moreover, the two methods may be extended to enable the differentiation of higher-order functions whose domains too are functions. In each case, these methods remedy a bug that would otherwise lead to the generation of incorrect results. An introduction to automatic differentiation techniques and to this bug are detailed below. 
     For expository purposes, we present this bug in the context of forward AD. However, the underlying issue can also manifest itself with other AD modes, including reverse AD of higher-order functions. The bug is insidious in that it can lead to production of incorrect results without warning. 
     Let   denote the true mathematical derivative operator. D is classically defined for first-order functions  →  in terms of limits, and thus this classical definition does not lend itself to direct implementation. 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             
                               𝔻 
                               ⁢ 
                               f 
                             
                             = 
                             
                               f 
                               ′ 
                             
                           
                         
                         
                           
                             where 
                             ⁢ 
                                 
                             
                               f 
                               ′ 
                             
                           
                         
                       
                     
                     ⁢ 
                     
                       ( 
                       x 
                       ) 
                     
                   
                   = 
                   
                     
                       lim 
                       
                         e 
                         → 
                         0 
                       
                     
                     
                       
                         
                           f 
                           ⁡ 
                           ( 
                           
                             x 
                             + 
                             ϵ 
                           
                           ) 
                         
                         - 
                         
                           f 
                           ⁡ 
                           ( 
                           x 
                           ) 
                         
                       
                       ϵ 
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     We seek to materialize   as a program construct  . We can view this classical limit definition as a specification of   and proceed to develop an implementation of  . Below, we use = to denote mathematical equality,   to denote definition of program constructs, and ⇒ to denote evaluation. 
     One can extend   to functions  →α, where: 
       α::= |α 1 →α 2   (2)
 
     Since by equation (2) any type α must be of the form α 1 → . . . →α n → , functions  →α 2  can be viewed as multivariate  →α 2 → . . . →α n →  whose first argument domain is   and whose range is  . We take   ƒ where ƒ: →α 2 → . . . →α n →  to be the partial derivative with respect to the first argument. 
     
       
         
           
             
               
                 
                   
                     𝔻 
                     ⁢ 
                     f 
                   
                   = 
                   
                     
                       ∂ 
                       
                         f 
                         ⁡ 
                         ( 
                         
                           
                             x 
                             1 
                           
                           , 
                           
                             x 
                             2 
                           
                           , 
                           … 
                               
                           , 
                           
                             x 
                             n 
                           
                         
                         ) 
                       
                     
                     
                       ∂ 
                       
                         x 
                         
                             
                           1 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     We will see below that a   can be implemented that appears to coincide with the specification   in equation (1) for functions  → , but then fails to coincide with the specification   in equation (3) for functions  →α. 
     Forward AD can be formulated as differential algebra. Its essence is as follows. 
     Complex numbers can be represented as pairs of real numbers that form an algebra over two-term polynomials a+bi where i 2 =−1. Arithmetic proceeds by simple rules, derived algebraically. 
       ( a+bi )+( c+di )=( a+c )+( b+d ) i   (4a)
 
       ( a+bi )( c+di )= ac +( ad+bc ) i+bdi   2 =( ac−bd )+( ad+bc ) i   (4b)
 
     Complex numbers can be implemented in a computer as ordered pairs (a, b), sometimes called Argand pairs. Since arithmetic over complex numbers is defined in terms of arithmetic over the reals, the above rules imply that computation over complex numbers is closed. 
     Dual numbers can be represented in the form a+b∈. In a dual number, the coefficient of e is called a perturbation or a tangent. These can similarly be viewed as an algebra over two-term polynomials where ∈ 2 =0 but ∈≠0. Arithmetic over dual numbers is again defined by simple rules derived algebraically. 
       ( a+b ∈)+( c+d ∈)=( a+c )+( b+d )∈  (5a)
 
       ( a+b ∈)( c+d ∈)= ac +( ad+bc )∈+ bd∈   2   =ac ( ad+bc )∈  (5b)
 
     Like complex numbers, dual numbers can be implemented in a computer as ordered pairs (a,b). Likewise, since arithmetic over dual numbers is defined in terms of arithmetic over the reals, the above rules imply that computation over dual numbers is closed. 
     The essence of forward AD is viewing dual numbers as truncated two-term power series. Since, following Taylor, ƒ(x 0 +x 1 ∈+O(∈ 2 ))=ƒ(x 0 )+x 1 ƒ′(x 0 )∈+O(∈ 2 ), applying ƒ to a dual number a+1∈ will yield a dual number ƒ(a)+ƒ′(a)∈. This leads to the following method for computing derivatives of functions ƒ: →  expressed as computer programs. 
     Particularly, given a programming language that supports dual numbers and arithmetic thereupon, ƒ′ at a point a can be computed by (Step 1) forming a +1∈, (Step 2) applying ƒ to a+1∈ to obtain a result ƒ(a)+ƒ′(a)∈, and (Step 3) extracting the tangent, ƒ′(a), from the result. 
     Step 2 constitutes a nonstandard interpretation of the arithmetic basis functions with equations (5a) and (5b). This can be implemented in various ways, for example, overloading or source-code transformation. Further, dual numbers can be represented in various ways, for example, as unboxed flattened values or as boxed values referenced through pointers. These different implementation strategies do not concern us here. While different implementation strategies have different costs, what we discuss applies to all strategies. 
     It is convenient to encapsulate Steps 1-3 as a higher-order function  :ƒ ƒ′. Indeed, that is one of the original motivations for the development of the lambda calculus. We can do this with the following code that implements  . 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           tg 
                           ⁢ 
                              
                           a 
                         
                         
                           = 
                           △ 
                         
                         0 
                       
                     
                     
                       
                                                                                    
                         
                           a 
                           : 
                           ℝ 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     6 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     tg 
                     ⁡ 
                     ( 
                     
                       a 
                       + 
                       
                         b 
                         ⁢ 
                         ϵ 
                       
                     
                     ) 
                   
                   
                     = 
                     △ 
                   
                   b 
                 
               
               
                 
                   ( 
                   
                     6 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     fx 
                   
                   
                     = 
                     △ 
                   
                   
                     tg 
                     ⁡ 
                     ( 
                     
                       f 
                       ⁡ 
                       ( 
                       
                         x 
                         + 
                         
                           1 
                           ⁢ 
                           ϵ 
                         
                       
                       ) 
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   
                     6 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     Here, x+1∈ denotes Step 1 above, that is, constructing a dual number, and tg (a+b∈) denotes Step 3 above, that is, extracting the tangent of a dual number. Equation (6a) handles the case where the output of ƒ is independent of the input x. 
     Forward AD provides certain complexity guarantees. Steps 1 and 3 take unit time. Step 2 introduces no more than a constant factor increase in both the space and time complexity of executing ƒ under a nonstandard interpretation. Thus, computing ƒ x and  ƒ x have the same space and time complexity. 
     In one implementation, dual numbers are tagged to avoid perturbation confusion. It is natural to nest application of  . Doing so would allow taking higher-order derivatives and, more generally, derivatives of functions that take derivatives of other functions. 
         (λ x  . . .  (λ y  . . . ) . . . ) . . .  (7)
 
     This can lead to perturbation confusion, yielding an incorrect result. The essence of perturbation confusion is that each invocation, calculation, and/or execution of   must perform its computation over a distinct differential algebra. While it is possible to reject programs that would exhibit perturbation confusion using static typing, and static typing can be used to yield the desired correct result in some cases with some user annotation, no static method is known that can yield the desired correct result in all cases without any annotation. It is possible, however, to get the correct result in all cases (except, as we shall see, when taking derivatives of functions whose ranges are functions) without user annotation, by redefining tg and   to tag dual numbers with distinct ∈s to obtain distinct differential algebras (or equivalently, distinct generators in a differential algebra) introduced by different invocations of  . We will indicate different tags by different subscripts on ∈, and use ε to denote a variable that is bound to an ∈. 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           tg 
                           ⁢ 
                              
                           ε 
                           ⁢ 
                              
                           a 
                         
                         
                           = 
                           △ 
                         
                         0 
                       
                     
                     
                       
                                                                                               
                         
                           a 
                           : 
                           ℝ 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     8 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     tg 
                     ⁢ 
                        
                     ε 
                     ⁢ 
                        
                     
                       ( 
                       
                         a 
                         + 
                         
                           b 
                           ⁢ 
                           ε 
                         
                       
                       ) 
                     
                   
                   
                     = 
                     △ 
                   
                   b 
                 
               
               
                 
                   ( 
                   
                     8 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       
                         
                           tg 
                           ⁢ 
                              
                           
                             ε 
                             1 
                           
                           ⁢ 
                              
                           
                             ( 
                             
                               a 
                               + 
                               
                                 b 
                                 ⁢ 
                                 
                                   ε 
                                   2 
                                 
                               
                             
                             ) 
                           
                         
                         
                           = 
                           △ 
                         
                         
                           
                             ( 
                             
                               tg 
                               ⁢ 
                                  
                               
                                 ε 
                                 1 
                               
                               ⁢ 
                                  
                               a 
                             
                             ) 
                           
                           + 
                           
                             
                               ( 
                               
                                 tg 
                                 ⁢ 
                                    
                                 
                                   ε 
                                   1 
                                 
                                 ⁢ 
                                    
                                 b 
                               
                               ) 
                             
                             ⁢ 
                             
                               ε 
                               2 
                             
                           
                         
                       
                     
                     
                       
                                    
                         
                           
                             ε 
                             1 
                           
                           ≠ 
                           
                             ε 
                             2 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     8 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     fx 
                   
                   
                     = 
                     △ 
                   
                   
                     fresh 
                     ⁢ 
                         
                     ε 
                     ⁢ 
                         
                     in 
                     ⁢ 
                         
                     tg 
                     ⁢ 
                        
                     
                       ε 
                       ⁡ 
                       ( 
                       
                         f 
                         ⁡ 
                         ( 
                         
                           x 
                           + 
                           
                             1 
                             ⁢ 
                             ε 
                           
                         
                         ) 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   
                     8 
                     ⁢ 
                     d 
                   
                   ) 
                 
               
             
           
         
       
     
     These redefine equations (6a), (6b), and (6c). Here, the tags are generated dynamically. Prior to this change, that is with only a single e, the values a and b in a dual number a+b∈ would be real numbers. With this change, that is with multiple ∈s, the values a and b in a dual number a+b∈ 1  can be dual numbers over ∈ 2  where ∈ 2 ∈ 1 . Such a tree of dual numbers will contain real numbers in its leaves and will contain a given ∈ only once along each path from the root to the leaves. Equation (8c) provides the ability to extract the tangent of an ∈ that might not be at the root of the tree. 
     This can be extended to higher-order functions whose range is a function, or in other words, to cases in which ƒ x returns a function g which takes an argument y. If one applies   to a function ƒ whose range is a function, then ƒ(x+1∈) in equation (8d) will yield a function  g  which takes the argument y. Thus, an invocation of  ƒ x yields a function  g  which takes the argument y and which itself performs a derivative calculation with respect to function ƒ x when invoked. It will not be possible to extract the tangent of this with tg as implemented by equations (8a), (8b), and (8c). The definition of tg can be augmented to handle this case by post-composition. 
         tgε g     ( tg ε)∘   g   g    is a function  (8e)
 
     However, this extension (alone) is flawed, as we proceed to demonstrate. Particularly, consider the following commonly occurring mathematical situation. We define an offset operator: 
         x : →( → )→ → 
 
         suf x   ƒ( x+u )  (9)
 
     The derivative of s at zero should be the same as the derivative operator, that is,  s 0= , since:    
     
       
         
           
             
               
                 
                   
                     
                       ( 
                       
                         ∀ 
                         f 
                       
                       ) 
                     
                     ⁢ 
                     
                       ( 
                       
                         ∀ 
                         y 
                       
                       ) 
                     
                     ⁢ 
                        
                     𝔻 
                     ⁢ 
                     s 
                     ⁢ 
                        
                     0 
                     ⁢ 
                     f 
                     ⁢ 
                     y 
                   
                   = 
                   
                     
                       
                         ∂ 
                         
                           ∂ 
                           u 
                         
                       
                       
                         
                           [ 
                           sufy 
                           ] 
                         
                         
                           u 
                           = 
                           0 
                         
                       
                     
                     = 
                     
                       
                         
                           
                             ∂ 
                             
                               ∂ 
                               u 
                             
                           
                           [ 
                           
                             f 
                             ⁡ 
                             ( 
                             
                               y 
                               + 
                               u 
                             
                             ) 
                           
                           ] 
                         
                         
                           u 
                           = 
                           0 
                         
                       
                       = 
                       
                         
                           
                             f 
                             ′ 
                           
                           ( 
                           y 
                           ) 
                         
                         = 
                         
                           𝔻 
                           ⁢ 
                           f 
                           ⁢ 
                           y 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 ⟺ 
               
               
                 
                                   
                   
                     { 
                     eta 
                     } 
                   
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       ( 
                       
                         ∀ 
                         f 
                       
                       ) 
                     
                     ⁢ 
                        
                     𝔻 
                     ⁢ 
                     s 
                     ⁢ 
                        
                     0 
                     ⁢ 
                     f 
                   
                   = 
                   
                     𝔻 
                     ⁢ 
                     f 
                   
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 ⟺ 
               
               
                 
                                   
                   
                     { 
                     eta 
                     } 
                   
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     𝔻 
                     ⁢ 
                     s 
                     ⁢ 
                        
                     0 
                   
                   = 
                   𝔻 
                 
               
               
                 
                   ( 
                   
                     10 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     Thus, if we define 
         ≙   s 0  (11)
 
     we would hope that  = . However, we exhibit an example where it does not. 
     We can compute  ( h)y for h: →  with simple reduction steps. Particularly,   is computed according to the steps: 
     
       
         
         
           
             {by(11)} 
           
         
       
    
           s 0  (12a)
         {by (8d)}       

       fresh ε in  tg ε( s (0+1ε))  (12b)
         {allocate a fresh tag ∈ 0 ; this is problematic; see discussion below}       

         tg∈   0 ( s (0+1∈ 0 ))  (12c)
         {by (9)}       

         tg∈   0 (λƒ·λ x ·(ƒ( x+ 1∈ 0 )))  (12d)
         {by (8e)}       

       ( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))  (12e)
         {postcompose}       

       λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 ))  (12f)
 
     Next,  ( h)y is computed according to the steps: 
       ( h)y
         {substitute (12f) for  }       
       (Δƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 )))((λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 ))) h ) y   (12g)
         {beta reduce}       

       (λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 )))(λ x·tg∈   0 ( h ( x+ 1∈ 0 ))) y   (12h)
         {beta reduce}       

       (λ x·tg∈   0 ((λ x·tg∈   0 ( h ( x+ 1∈ 0 )))( x+ 1∈ 0 ))) y   (12i)
         {beta reduce}       

         tg∈   0 ((λ x·tg∈   0 ( h ( x+ 1∈ 0 )))( y+ 1∈ 0 ))  (12j)
         {beta reduce}       

         tg∈   0 ( tg∈   0 ( h (( y+ 1∈ 0 )+1∈ 0 )))  (12k)
         {add dual numbers}       

         tg∈   0 ( tg∈   0 ( h ( y+ 2∈ 0 )))  (12l)
         {apply h to a dual number}       

         tg∈   0 ( tg∈   0 ( h ( y )+2 h ′( y )∈ 0 ))  (12m)
         {by (8b)}       

         tg∈   0 (2 h ′( y ))  (12n)
         {by (8a)}       

       0  (12o)
 
     As can be seen, this went wrong, yielding 0 instead of h″(y). Particularly: 
         (   h ) y⇒ 0≠ (   h ) y=h ″( y )  (13)
 
     The process of allocating a fresh tag in step (12d) was problematic. The proper way to handle such fresh tag allocation might be to use nominal logic, perhaps in a dependent-type-theoretic variant. Below, we offer alternate mechanisms that are suitable for use in programming-language implementations that lack type systems that support first class names and binding. 
     This is not an artificial example. It is quite natural to construct an x-axis differential operator and apply it to a two-dimensional function twice, along the x and then y axis directions, by applying the operator, flipping the axes, and applying the operator again, thus creating precisely this sort of cascaded use of a defined differential operator. 
     This incorrect result was due to the tag ∈ 0  being generated exactly once, in equation (12b), when   was calculated from  s 0 as equations (12a)-(12f) using the definition of equation (11). The invocation  s 0 is the point at which a fresh tag is introduced; early instantiation can result in reuse of the same tag in logically distinct derivative calculations. Here, the first derivative and the second derivative become confused at equation (12l). We have two nested applications of tg for ∈ 0 , but for correctness these should be distinctly tagged: ∈ 0  versus ∈ 1 . 
     This can be accomplished by making two copies of   by evaluating   s 0 twice. Performing an analogous computation with two copies of   yields the correct result. Particularly,    0  is computed according to the steps: 
     
       
       0  
         
         
           
             {repeat (12a)} 
           
         
       
    
           s 0  (14a)
         {repeat (12b)}       

       fresh ε in  tg ε( s (0+1ε))  (14b)
         {repeat (12c)}       

         tg∈   0 ( s (0+1∈ 0 ))  (14c)
         {repeat (12d)}       

         tg∈   0 (λƒ·λ x ·(ƒ( x+ 1∈ 0 )))  (14d)
         {repeat (12e)}       

       ( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))  (14e)
         {repeat (12f)}       

       λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 ))  (14f)
 
     Next,    1  is computed according to the steps: 
     
       
       1  
         
         
           
             {repeat (12a)} 
           
         
       
    
           s 0  (14g)
         {repeat (12b)}       

       fresh ε in  tg ε( s (0+1ε))  (14h)
         {repeat (12c)}       

         tg∈   1 ( s (0+1∈ 1 ))  (14i)
         {repeat (12d)}       

         tg∈   1 (λƒ·λ x ·(ƒ( x+ 1∈ 1 )))  (14j)
         {repeat (12e)}       

       ( tg∈   1 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 1 )))  (14k)
         {repeat (12f)}       

       λƒ·λ x·tg∈   1 (ƒ( x+ 1∈ 1 ))  (14l)
 
     Finally,    0 (   1 h)y is computed according to the steps: 
         0 (   1 h)y
         {substitute (14f) and (14l) for  }       
       (λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 )))((λƒ·λ x·tg∈   1 (ƒ( x+ 1∈ 1 ))) h ) y   (14m)
         {beta reduce}       

       (λƒ·λ x·tg∈   0 (ƒ( x+ 1∈ 0 )))(λ x·tg∈   1 ( h ( x+ 1∈ 1 ))) y   (14n)
         {beta reduce}       

       (λ x·tg∈   0 ((λ x·tg∈   1 ( h ( x+ 1∈ 1 )))( x+ 1∈ 0 ))) y   (14o)
         {beta reduce}       

         tg∈   0 ((λ x·tg∈   1 ( h ( x+ 1∈ 1 )))( y+ 1∈ 0 ))  (14p)
         {beta reduce}       

         tg∈   0 ( tg∈   1 ( h (( y+ 1∈ 0 )+1∈ 1 )))  (14q)
         {apply h to a dual number}       

         tg∈   0 ( tg∈   1 ( h ( y+ 1∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))  (14r)
         {apply h to a dual number}       

         tg∈   0 ( tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))  (14s)
         {apply h to a dual number}       

         tg∈   0 ( tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+( h ′( y )+ h ″( y )∈ 0 )∈ 1 ))  (14t)
         {by (8b)}       

         tg∈   0 ( h ′( y )+ h ″( y )∈ 0 )  (14u)
 
     {by (8b)} 
         h ″( y )  (14v)
 
     Here, equation (14r) corrects the mistake in equation (12l). 
     However, this is tantamount to requiring the user to manually write 
       let    0       s 0 
       in let    1       s 0 
       in    0 (   1   h ) y   (15)
 
       instead of: 
       let        s 0 
       in  (   h ) y   (16)
 
     This should not be necessary since if   correctly implemented  . Particularly,    0  and    1  should be equivalent when   correctly implements  . 
     The essence of the bug is that the implementation of   in equation (8d) only generates a distinct ∈ for each invocation   ƒ x, but a distinct ∈ is instead needed for each derivative calculation. In the first-order case, when ƒ: → , these are equivalent. Each invocation   ƒ x leads to a single derivative calculation. But in the higher-order case, when ƒ returns a function g, an invocation   ƒ x yields  g  which performs a derivative calculation when invoked. Since  g  can be invoked multiple times, each such invocation will perform a distinct derivative calculation and needs a distinct ε. 
     Methods for more accurately implementing   are described below. The methods described below advantageously extend AD to support functions whose domains and/or ranges are functions and, importantly, remedy the insidious bug described above. Thus, these methods enable AD that is completely general and which can be applied in an unrestricted fashion to correctly compute the derivative of all programs that compute differentiable mathematical functions. This includes application to functions whose domain and/or ranges include the entire space of data types supported by programming languages, including not only aggregates but also functions. 
     First Method of Improved Automatic Differentiation of Higher-Order Functions 
     A first method for resolving the bug described above is to eta expand the definition of  . Such eta expansion would need to be conditional on the return type of ƒ. 
     
       
         
           
             
               
                 
                   
                     
                       1 
                     
                     : 
                     
                       
                         ( 
                         
                           ℝ 
                           → 
                           ℝ 
                         
                         ) 
                       
                       → 
                       ℝ 
                       → 
                       ℝ 
                     
                   
                   ⁢ 
                     
                   
 
                   
                     
                       
                         1 
                       
                       f 
                       ⁢ 
                       
                         x 
                         1 
                       
                     
                     
                       = 
                       △ 
                     
                     
                       fresh 
                       ⁢ 
                           
                       ε 
                       ⁢ 
                           
                       in 
                       ⁢ 
                           
                       tg 
                       ⁢ 
                          
                       
                         ε 
                         ⁡ 
                         ( 
                         
                           f 
                           ⁡ 
                           ( 
                           
                             
                               x 
                               1 
                             
                             + 
                             
                               1 
                               ⁢ 
                               ε 
                             
                           
                           ) 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     17 
                     ⁢ 
                     a 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       2 
                     
                     : 
                     
                       
                         ( 
                         
                           ℝ 
                           → 
                           
                             α 
                             2 
                           
                           → 
                           ℝ 
                         
                         ) 
                       
                       → 
                       ℝ 
                       → 
                       
                         α 
                         2 
                       
                       → 
                       ℝ 
                     
                   
                   ⁢ 
                     
                   
 
                   
                     
                       
                         2 
                       
                       f 
                       ⁢ 
                       
                         x 
                         1 
                       
                       ⁢ 
                       
                         x 
                         2 
                       
                     
                     
                       = 
                       △ 
                     
                     
                       fresh 
                       ⁢ 
                           
                       ε 
                       ⁢ 
                           
                       in 
                       ⁢ 
                           
                       tg 
                       ⁢ 
                          
                       
                         ε 
                         ⁡ 
                         ( 
                         
                           
                             f 
                             ⁡ 
                             ( 
                             
                               
                                 x 
                                 1 
                               
                               + 
                               
                                 1 
                                 ⁢ 
                                 ε 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                             x 
                             2 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     17 
                     ⁢ 
                     b 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               
                 
                   
                     
                       3 
                     
                     : 
                     
                       
                         ( 
                         
                           ℝ 
                           → 
                           
                             α 
                             2 
                           
                           → 
                           
                             α 
                             3 
                           
                           → 
                           ℝ 
                         
                         ) 
                       
                       → 
                       ℝ 
                       → 
                       
                         α 
                         2 
                       
                       → 
                       
                         α 
                         3 
                       
                       → 
                       ℝ 
                     
                   
                   ⁢ 
                     
                   
 
                   
                     
                       
                         3 
                       
                       f 
                       ⁢ 
                       
                         x 
                         1 
                       
                       ⁢ 
                       
                         x 
                         2 
                       
                       ⁢ 
                       
                         x 
                         3 
                       
                     
                     
                       = 
                       △ 
                     
                     
                       fresh 
                       ⁢ 
                           
                       ε 
                       ⁢ 
                           
                       in 
                       ⁢ 
                           
                       tg 
                       ⁢ 
                          
                       
                         ε 
                         ⁡ 
                         ( 
                         
                           
                             f 
                             ⁡ 
                             ( 
                             
                               
                                 x 
                                 1 
                               
                               + 
                               
                                 1 
                                 ⁢ 
                                 ε 
                               
                             
                             ) 
                           
                           ⁢ 
                           
                             x 
                             2 
                           
                           ⁢ 
                           
                             x 
                             3 
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                     
                   ⋮ 
                 
               
               
                 
                   ( 
                   
                     17 
                     ⁢ 
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
     It should be appreciated the “eta expansion” refers to a process of syntactically expanding an expression. In particular, the eta expansion of an expression E refers to the transformation of the expression into the expanded expression λv·Ev, where the variable v is fresh and E:α→β. The expression E must be a function. Eta expansion delays the evaluation of E until the function is applied, and will re-evaluate E each time the function is applied. 
     As noted above, if one applies   to a higher-order function ƒ x that returns a function g which takes an argument y, then ƒ(x+1∈) in equation (8d) will yield a function  g  which takes the argument y. Thus, an invocation of   ƒ x yields a function g which takes the argument y and which itself performs a derivative calculation with respect to function ƒ x when invoked. It will not be possible to extract the tangent of this with tg as implemented by equations (8a), (8b), and (8c). However, with such eta expansion conditioned on the return type of ƒ, this situation can be avoided and equation (8e) is not needed, because the appropriate variant of   should only be invoked in a context that contains all arguments necessary to subsequently allow the call to tg in that invocation of   to yield to a non-function-containing value. This seemingly infinite set of    i  and associated definitions can be formulated as a single   with polymorphic recursion. 
         ƒ x     λy ·( (λ x ·(ƒ xy )) x ) (ƒ x ) is a function  (18a)
 
         ƒ x   fresh ε in  tg ε(ƒ( x+ 1ε)) (ƒ x ) is not a function  (18b)
 
     We can see that this resolves the bug in equations (12a)-(12o) and accomplishes the desiderata in equations (14a)-(14l) without making two copies of  . Particularly,   is computed according to the steps: 
     
       
         
         
           
             {by (11)} 
           
         
       
    
           s 0  (19a)
         {by (18a)}       

       λ y ·( (λ x ·( sxy ))0)  (19b)
 
     Next,  ( h)y is computed according to the steps: 
       ( h)y
         {substitute (19b) for  }       
       (λ y ·( (λ x ·( sxy ))0))((λ y ·( (λ x ·( sxy ))0)) h ) y   (19c)
         {beta reduce}       

       (λ y ·( (λ x ·( sxy ))0))( (λ x ·( sxh ))0) y   (19d)
         {beta reduce}       

       ( (λ x ·( sx ( (λ x·sxh )0)))0) y   (19e)
         {by (8d)}       

       (fresh ε in  tg ε((λ x ·( sx ( (λ x ·( sxh ))0)))(0+1ε))) y   (19f)
         {allocate a fresh tag ∈ 0 }       

       ( tg∈   0 ((λ x ·( sx ( (λ x ·( sxh ))0)))(0+1∈ 0 ))) y   (19g)
         {beta reduce}       

       ( tg∈   0 ( s (0+∈ 0 )( (λ x ·( sxh ))0))) y   (19h)
         {by (8d)}       

       ( tg∈   0 ( s (0+1∈ 0 )(fresh ε in  tg ε((λ x ·( sxh ))(0+1ε))))) y   (19i)
         {allocate a fresh tag ∈ 1 }       

       ( tg∈   0 ( s (0+1∈ 0 )( tg∈   1 (λ x ·( sxh ))(0+1∈ 1 ))))) y   (19j)
         {beta reduce}       

       ( tg∈   0 ( s (0+1∈ 0 )( tg∈   1 ( s (0+1∈ 1 ) h )))) y   (19k)
         {by (9)}       

       ( tg∈   0 ( s (0+1∈ 0 )( tg∈   1 (λ x ·( h ( x +(0+1∈ 1 ))))))) y   (19l)
         {by (8e)}       

       ( tg∈   0 ( s (0+1∈ 0 )( tg∈   1 )∘(λ x ·( h ( x +(0+1∈ 1 )))))) y   (19m)
         {postcompose}       

       ( tg∈   0 ( s (0+1∈ 0 )(λ x ·( tg∈   1 ( h ( x +(0+1∈ 1 ))))))) y   (19n)
         {by (9)}       

       ( tg∈   0 (λ x ·((λ x ·( tg∈   1 ( h ( x +(0+1∈ 1 )))))( x +(0+1∈ 0 ))))) y   (19o)
         {beta reduce}       

       ( tg∈   0 (λ x ·( tg∈   1 ( h (( x +(0+1∈ 0 ))+(0+1∈ 1 ))))) y   (19p)
         {by (8e)}       

       ( tg∈   0 )∘(λ x ·( tg∈   1 ( h (( x +(0+1∈ 0 ))+(0+1∈ 1 ))))) y   (19q)
         {postcompose}       

       (λ x ·( tg∈   0 ( tg∈   1 ( h (( x +(0+1∈ 0 ))+(0+1∈ 1 )))))) y   (19r)
         {beta reduce}       

         tg∈   0 ( tg∈   1 ( h (( y +(0+1∈ 0 ))+(0+1∈ 1 ))))  (19s)
         {add dual numbers}       

         tg∈   0 ( tg∈   1 ( h (( y+ 1∈ 0 )+(0+1∈ 1 )))))  (19t)
         {add dual numbers}       

         tg∈   0 ( tg∈   1 ( h (( y+ 1∈ 0 ))+1∈ 1 )))  (19u)
         {same as (14r)}       

         tg∈   0 ( tg∈   1 ( h ( y+ 1∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))  (19v)
         {same as (14s)}       

         tg∈   0 ( tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))  (19w)
         {same as (14t)}       

         tg∈   0 ( tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+( h ′( y )+ h ″( y )∈ 0 )∈ 1 ))  (19x)
         {same as (14u)}       

         tg∈   0 ( h ′( y )+ h ″( y )∈ 0 )  (19y)
         {same as (14v)}       

         h ″( y )  (19z)
 
     Here, the allocation of a fresh tag is delayed from equation (19b) and is performed twice, in equations (19g) and (19j), allowing equation (19v) to correct the mistake in equation (12l), just like equation (14r). 
     It should be appreciated that this disclosure only considers a space of types that includes scalar reals and functions but not aggregates (exclusive of dual numbers). Complications can arise when extending the space of types to include aggregates. The example below illustrates that the above mechanism works with functions that return Church-encoded aggregates. 
       ( a,d ) m     m a d   (20a)
 
         fst c     c (λ a ·(λ d·a )  (20b)
 
         snd c     c (λ a ·(λ d·d ))  (20c)
 
         tu   ( e   u×u ,(λƒ·(λ x ·(ƒ x+u ))))  (20d)
 
           t 1⇒ t ′(1)  (20e)
 
         p       t 0  (20f)
 
         fst p⇒ 0  (20g)
 
             snd p   (20h)
 
         (  exp)1⇒ e   (20i)
 
     With a function that returned native aggregates, one would need to emulate the behavior that occurs with Church-encoded aggregates on native aggregates by delaying derivative calculation, with the associated tag allocation and tg applied to the native retuned aggregate, until an accessor is applied to that aggregate. Consider  t 0 where t: →( ×(( → )→ )) as above. One could not perform the derivative calculation when computing the value p returned by  t 0. One would have to delay until applying an accessor to p. If one accessed the first element of p, one would perform the derivative calculation, with the associated tag allocation, at the time of access. But if one accessed the second element of p, one would have to further delay the derivative calculation, with the associated tag allocation, until that second element was invoked. This could require different amounts of delay that might be incompatible with some static type systems. 
     Moreover, it will be appreciated that a type system or other static analysis mechanism that is unable to handle the unbounded polymorphism of equations (17a), (17b), (17c), . . . or infer the “is [not] a function” side conditions of equations (18a) and (18b), achieving completeness might require run-time evaluation of the side conditions. This could involve calling ƒ twice, once to determine its return type and once to do the eta-expanded derivative calculation, and lead to exponential increase in asymptotic time complexity. 
     Finally, the solution may break sharing in curried functions, even with a type system or other static analysis mechanism that is able to eliminate the run-time evaluation of “is [not] a function” side conditions. Consider 
         gx   let  t     ƒx  in λ p·pt   (21)
 
       invoked in: 
         hx   let  c     gx  in ( c (λ t·t ))+( c (λ t ·(λ u·t×u ))π)  (22)
 
     The programmer would expect h 8 to call ƒ once in the calculation of the temporary t=ƒ 8. And indeed this is what would occur in practice. Now consider   h 8. The strategy discussed above would (in the absence of memoization or similar heroic measures) end up calculating ƒ 8 twice, as the delayed tag allocation would end up splitting into two independent tag allocations with each independently redoing the calculation. This violates the constant-factor-overhead complexity guarantee of forward AD, imposing, in the worst case, exponential overhead. 
     Second Method of Improved Automatic Differentiation of Higher-Order Functions 
     A second method for resolving the bug described above is to wrap  g  with tag substitution to guard against tag collision. Particularly, as noted above, if one applies   to a higher-order function ƒ x that returns a function g which takes an argument y, then ƒ(x+1ε) in equation (8d) will yield a function  g  which takes the argument y. Thus, an invocation of  ƒ x yields a function  g  which takes the argument y and which itself performs a derivative calculation with respect to function ƒ x when invoked. It will not be possible to extract the tangent of this with tg as implemented by equations (8a), (8b), and (8c). The definition of tg can be augmented to handle this and to also resolve the bug in equation (8e) by replacing equation (8e) with: 
         tgε   1     g y   fresh ε in ([ε 1 /ε]∘( tgε   1 )∘   g   ∘[ε/ε 1 ]) y  g    is  a  function  (23)
 
     Here [ε 1 /ε 2 ] x substitutes ε 1  for ε 2  in x. In a language with opaque closures, tag substitution must operate on functions by appropriate pre- and post-composition. 
       [ε 1 /ε 2 ] a     a a:         (24a)
 
       [ε 1 /ε 2 ]( a+bε   2 )   a+bε   1   (24b)
 
       [ε 1 /ε 2 ]( a+b ε) ([ε 1 /ε 2 ] a )+([ε 1 /ε 2 ] b )ε ε≠ε 2   (24c)
 
       [ε 1 /ε 2 ]   g y   fresh ε in ([ε 2 /ε]∘[ε 1 /ε 2 ]∘   g   ∘[ε/ε 2 ]) y  g    is  a  function  (24d)
 
     The intent of equation (24d) is to substitute ε 1  for ε 2  in values closed-over in  g . An ε 2  in the output of  g  can result either from closed-over values and/or input values. We want to substitute for instances of ε 2  in the output that result from the former but not the latter. This is accomplished by substituting a fresh tag for instances of ε 2  in the input and substituting them back at the output to preserve the extensional behavior of  g . Equation (23) operates in a similar fashion. The intent of equation (23) is to extract the coefficient of instances of ε 1  in the output of  g  that result from closed-over values, not input values. This is accomplished by substituting a fresh tag for instances of ε 1  in the input and substituting them back at the output to preserve the extensional behavior of  g . 
     We can see that this also resolves the bug in equations (12a)-(12o) and accomplishes the desiderata in equations (14a)-(14l) without making two copies of  . Particularly,   is computed according to the steps: 
     
       
         
         
           
             {by (11)} 
           
         
       
    
           s 0  (25a)
         {by (8d)}       

       fresh ε in  tg ε( s (0+1ε))  (25b)
         {allocate a fresh tag ∈ 0 }       

         tg∈   0 ( s (0+1∈ 0 ))  (25c)
         {by (9)}       

         tg∈   0 (λƒ·λ x ·(ƒ( x+ 1∈ 0 )))  (25d)
         {by (23)}       

       Δ y ·(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) y )  (25e)
 
     Next,  ( h)y is computed according to the steps: 
       ( h)y
         {substitute (25e) for  }       
       λ y ·(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) y )(λ y ·(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) y ) h ) y   (25f)
         {beta reduce}       

       λ y ·(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) y )(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ(+1∈ 0 )))∘[ε/∈ 0 ]) h ) y   (25g)
         {beta reduce}       

       (fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ])(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) h )) y   (25h)
         {allocate a fresh tag ∈ 1 }       

       (([∈ 0 /∈ 1 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[∈ 1 /∈ 0 ])(fresh ε in ([∈ 0 /ε]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[ε/∈ 0 ]) h )) y   (25i)
         {allocate a fresh tag ∈ 2 }       

       (([∈ 0 /∈ 1 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[∈ 1 /∈ 0 ])(([∈ 0 /∈ 2 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[∈ 2 /∈ 0 ]) h )) y   (25j)
         {substitute ∈ 2  for ∈ 0 , which leaves h unchanged since it can&#39;t close over the freshly allocated tags}       

       (([∈ 0 /∈ 1 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[∈ 1 /∈ 0 ])(([∈ 0 /∈ 2 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))) h )) y   (25k)
         {beta reduce and postcompose}       

       (([∈ 0 /∈ 1 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 )))∘[∈ 1 /∈ 0 ])(λ x ·([∈ 0 /∈ 2 ]( tg∈   0 ( h ( x+ 1∈ 0 )))))) y   (25l)
         {substitute ∈ 1  for ∈ 0 }       

       (([∈ 0 /∈ 1 ]∘( tg∈   0 )∘(λƒ·λ x ·(ƒ( x+ 1∈ 0 ))))(λ x ·([∈ 1 /∈ 2 ]( tg∈   1 ( h ( x+ 1∈ 1 )))))) y   (25m)
         {beta reduce and postcompose}       

       (λ x ·([∈ 0 /∈ 1 ]( tg∈   0 ((λ x ·([∈ 1 /∈ 2 ]( tg∈   1 ( h ( x+ 1∈ 1 )))))( x+ 1∈ 0 ))))) y   (25n)
         {beta reduce}       

       [∈ 0 /∈ 1 ]( tg∈   0 ((λ x ·([∈ 1 /∈ 2 ]( tg∈   1 ( h ( x+ 1∈ 1 )))))( y+ 1∈ 0 )))  (25o)
         {beta reduce}       

       [∈ 0 /∈ 1 ]( tg∈   0 ([∈ 1 /∈ 2 ]( tg∈   1 ( h (( y+ 1∈ 0 )+1∈ 1 )))))  (25p)
         {apply h to a dual number}       

       [∈ 0 /∈ 1 ]( tg∈   0 ([∈ 1 /∈ 2 ]( tg∈   1 ( h ( y+ 1∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))))  (25q)
         {apply h to a dual number}       

       [∈ 0 /∈ 1 ]( tg∈   0 ([∈ 1 /∈ 2 ] tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+ h ′( y+ 1∈ 0 )∈ 1 ))))   (25r)
         {apply h to a dual number}       

       [∈ 0 /∈ 1 ]( tg∈   0 ([∈ 1 /∈ 2 ] tg∈   1 (( h ( y )+ h ′( y )∈ 0 )+ h ′( y )+ h ″( y )∈ 0 )∈ 1 ))))  (25s)
         {by (8b)}       

       [∈ 0 /∈ 1 ]( tg∈   0 ([∈ 1 /∈ 2 ]( h ′( y )+ h ″( y )∈ 0 )))  (25t)
         {substitute ∈ 1  for ∈ 2 }       

       [∈ 0 /∈ 1 ]( tg∈   0 ( h ′( y )+ h ″( y )∈ 0 ))  (25u)
         {by (8b)}       

       [∈ 0 /∈ 1 ] h ″( y )  (25v)
         {substitute ∈ 0  for ∈ 1 }       

         h ″( y )  (25w)
 
     Equations (25k) and (25m) are abbreviated as they really use equation (24d). Here, the tag substitution in equation (25m) allows equation (25q) to correct the mistake in equation (12l), just like equation (14r). 
     It should be appreciated that this solution can present problems when implemented as user code in a pure language. In the presence of aggregates, unless care is taken, the computational burden of tag substitution can violate the complexity guarantees of forward AD. The call to tg in Step 3 might take longer than unit time as tag substitution must potentially traverse an aggregate of arbitrary size. When that aggregate shares substructure, a careless implementation might traverse such shared substructure multiple times, leading to potential exponential growth in time complexity. Moreover, a careless implementation might copy shared substructure multiple times, leading to potential exponential growth in space complexity. 
     In one embodiment, laziness, memoization, and hash-consing are utilized to help solve this, but it can be tricky to employ such in a fashion that preserves the requisite time and space complexity guarantees of forward AD, particularly in a pure or multithreaded context. Moreover, it will be appreciated that laziness, memoization, and hash-consing might not completely eliminate the problem. First, some languages like PYTHON and SCHEME lack the requisite pervasive default laziness. Failure to explicitly code the correct portions of user code as lazy in an eager language can break the complexity guarantees in subtle ways. But there are subtle issues even in languages like HASKELL with the requisite pervasive default laziness, and even when laziness is correctly introduced manually in eager languages. One is that memoization and hash-consing implicitly involve a notion of equality. But it is not clear what notion of equality to use, especially with “gensym” and potential alpha equivalence. One might need eq?, pointer or intentional equivalence, rather than equal?, structural or extensional equivalence, and all of the impurity that this introduces. Further, memoization and hash-consing might themselves be a source of a new kind of perturbation confusion if tags can persist. One would then need to substitute the memorized tags or the hash-cons cache. Beyond this, memoization and hash-consing could break space complexity guarantees unless the cache were flushed. It is not clear when/where to flush the cache, and even whether there is a consistent place to do so. There might be inconsistent competing concerns. Finally, many systems don&#39;t provide the requisite hooks to do all of this. One would need weak pointers and finalization. 
     The above difficulties only arise when implementing tag substitution as user code in a pure language. The opacity of closures necessitates implementing tag substitution on functions via pre- and post-composition (24d). The complexity guarantees of forward AD could be maintained if the substitution mechanism [ε 1 /ε 2 ]x were implemented so that it (a) did not traverse shared substructure multiple times, (b) copied shared sub structure during renaming in a fashion that preserved structure sharing, and (c) could apply to closures, by accessing, copying, renaming, and reclosing around the environments inside closures, without resorting to pre- and post-composition. This could be accomplished either by including the [ε 1 /ε 2 ]x mechanism as a primitive in the implementation, or by providing other lower-level primitives out of which it could be fashioned. One such mechanism is map-closure, the ability to reflectively access and modify closure environments. 
     Extension of the Methods to Higher-Order Functions Whose Ranges and Domains are Functions 
     The definition of equation (3) only extends  , and the mechanisms of the first method and second method described above only extend   to higher-order functions  →α whose ranges are functions. Differential geometry provides the framework for extending   to functions α 1 →α 2  whose domains too are functions. 
     Differential geometry concerns itself with differentiable mappings between manifolds, where intuitively a manifold is a surface along which points can move smoothly, like the surface of a sphere or the space of n×n rotation matrices. Given a point x, called a primal (value), on a manifold α, we can consider infinitesimal perturbations of x. The space of such perturbations is a vector space called a tangent space, denoted by T x α. This is a dependent type, dependent on the primal x. A particular perturbation, an element x′ of the tangent space, is called a tangent (value). A pair (x, x′) of a primal and tangent value is called a bundle (value), which are members of a bundle space Tα=Σ x:α {x}×T x α. Bundles generalize the notion of dual numbers. So, if x has type α, for some α, the tangent x′ has type T x α, and they can be bundled together as (x+x′∈) which has type Tα. 
     The machinery of differential geometry defines T x α for various manifolds and spaces α. For function spaces α→β, where ƒ is of type α→β, T ƒ (α→β)=(a:α)→T ƒ(a) β and T(α→β)=α→Tβ. The function bundle (x:α)(x′:T x α) (x:x′):Tα constructs a bundle from a primal and a tangent, and the function tangent (x,x′):Tα x′:T x α extracts a tangent from a bundle. Differential geometry provides a push forward operator that generalizes the notion of a univariate derivative from functions ƒ of type  →  to functions ƒ of type α→β. 
         pf :(α→β)→( Tα→T β)  (26)
 
     This augments the original mapping (a:α)→β to also linearly map a tangent T a α of the input a to a tangent T ƒ(a) β of the output ƒ(a). 
     Here we sketch how to materialize differential geometry as program constructs to generalize   to functions α 1 →α 2  whose domains (and ranges) are functions. We first note that: 
         ƒ  x =tangent( pf ƒ(bundle×1))  (27)
 
     This only applies when x:  because of the constant 1. We can generalize this to a directional derivative: 
         ƒ  x x ′=tangent( pf ƒ(bundle  x x ′))  (28)
 
     This further generalizes to x of any type. With this,   becomes a special case of  : 
         ƒ  x=     ƒ× 1  (29)
 
     To materialize   in equation (28), we need to materialize tangent, pf, and bundle. The definition of tg in equations (8a)-(8c) and equation (8e) materializes tangent with the first method that utilizes eta expansion. Likewise, definition of tg in equations (8a)-(8c) and equation (23) materializes tangent with the second method that utilizes tag substitution. The nonstandard interpretation of the arithmetic basis functions sketched in equations (5a) and (5b) materializes pf by lifting a computation on real numbers to a computation on dual numbers. All that remains is to materialize bundle. So far, we have been simply writing this as Step 2, a map from a to a+1∈ or a map from x to x+1ε in equation (8d). This only works for numbers, not functions. 
     With the framework of the first method that utilizes eta expansion, we can extend this to functions: 
       bun ε  x x′≙x+x′ε x  and  x ′ are not functions  (30a)
 
       bun εƒ ƒ′ y ≙bunε(ƒ  y )(ƒ′  y ) ƒ and ƒ′ are functions  (30b)
 
     The post-composition in equation (30b) is analogous to that in equation (8e). 
     With the framework of the second method that utilizes tag substitution, instead of to equation (30b), we use the alternative: 
       bunε 1   ƒ ƒ′ y ≙fresh ε ƒ and ƒ′ are functions
 
       in [ε 1 /ε](bunε 1 (ƒ([ε/ε 1 ] y ))(ƒ′([ε/ε 1 ] y )))  (31)
 
     The additional tag substitution in equation (31) is analogous to that in equation (23). 
     With this, we can now materialize j in the framework of the first method that utilizes eta expansion: 
         ƒ xx′     λy ,( (λ x ,(ƒ xy )) xx ′) (ƒ x ) is a function  (32a)
 
         ƒ xx ′ fresh ε in  tg ε(ƒ(bunε xx ′)) (ƒ x ) is not a function  (32b)
 
     The equations (32a) and (32b) are analogous to equations (18a) and (18b). 
     Likewise, we can now materialize   in the framework of the second that utilizes tag substitution: 
         ƒ xx ′ fresh ε in  tg ε(ƒ(bunε xx ′))  (33)
 
     The equation (33) is analogous to equation (8d). 
     With this,   becomes a special case of  : 
         ƒ x       ƒ× 1  (34)
 
     There is a crucial difference, however, between bundle and tangent and the corresponding materializations bun and tg. The former do not take e as an argument. This allows them to be used as distinct notational entities. In contrast, bun and tg must take the same e as an argument, this tag must be fresh, and it should not be used anywhere else. Thus, it should not escape, except in ways that are protected by tag substitution. This motivates creation of the   construct. There is no corresponding standard   construct in differential geometry; we created it just to describe the intended meaning of  . 
     Automatic Differentiation System 
       FIG. 1  shows a block diagram of an exemplary embodiment of an automatic differentiation system  100 . The automatic differentiation system  100  advantageously utilizes the methods described above to provide automatic differentiation (AD) that accurately supports functions whose domains and/or ranges are functions. These methods advantageously enable AD that is completely general and which can be applied in an unrestricted fashion to correctly compute the derivative of all programs that compute differentiable mathematical functions. This includes application to functions whose domain and/or ranges include the entire space of data types supported by programming languages, including not only aggregates but also functions. Moreover, as described in detail above, these methods advantageously remedy an insidious bug that would otherwise lead to incorrect results. 
     In the illustrated exemplary embodiment, the automatic differentiation system  100  comprises at least one processor  102 , at least one memory  104 , a communication module  106 , a display screen  108 , and a user interface  110 . However, it will be appreciated that the components of the automatic differentiation system  100  shown and described are merely exemplary and that the automatic differentiation system  100  may comprise any alternative configuration. Particularly, the automatic differentiation system  100  may comprise any computing device such as a desktop computer, a laptop, a smart phone, a tablet, or other personal electronic device. Thus, the automatic differentiation system  100  may comprise any hardware components conventionally included in such computing devices. 
     The memory  104  is configured to store data and program instructions that, when executed by the at least one processor  102 , enable the automatic differentiation system  100  to perform various operations described herein. The memory  104  may be of any type of device capable of storing information accessible by the at least one processor  102 , such as a memory card, ROM, RAM, hard drives, discs, flash memory, or any of various other computer-readable medium serving as data storage devices, as will be recognized by those of ordinary skill in the art. Additionally, it will be recognized by those of ordinary skill in the art that a “processor” includes any hardware system, hardware mechanism or hardware component that processes data, signals or other information. Thus, the at least one processor  102  may include a central processing unit, graphics processing units, multiple processing units, dedicated circuitry for achieving functionality, programmable logic, or other processing systems. Additionally, it will be appreciated that, although the automatic differentiation system  100  is illustrated as single device, the automatic differentiation system  100  may comprise several distinct processing systems  40  that work in concert to achieve the functionality described herein. 
     The communication module  106  may comprise one or more transceivers, modems, processors, memories, oscillators, antennas, or other hardware conventionally included in a communications module to enable communications with various other devices, such as the drone  30 . In at least some embodiments, the communication module  106  includes a Wi-Fi module configured to enable communication with a Wi-Fi network and/or Wi-Fi router (not shown). In further embodiments, the communications modules  46  may further include a Bluetooth® module, an Ethernet adapter and communications devices configured to communicate with wireless telephony networks. 
     The display screen  108  may comprise any of various known types of displays, such as LCD or OLED screens. In some embodiments, the display screen  108  may comprise a touch screen configured to receive touch inputs from a user. The user interface  110  may suitably include a variety of devices configured to enable local operation of the automatic differentiation system  100  by a user, such as a mouse, trackpad, or other pointing device, a keyboard or other keypad, speakers, and a microphone, as will be recognized by those of ordinary skill in the art. Alternatively, in some embodiments, a user may operate the automatic differentiation system  100  remotely from another computing device which is in communication therewith via the communication module  106  and has an analogous user interface. 
     The program instructions stored on the memory  104  include an automatic differentiation program  112 . As discussed in further detail below, the processor  102  is configured to execute the automatic differentiation program  112  to (i) receive first program code that defines a function, in particular a higher-order function whose range and/or domain are functions and (ii) determine or evaluate a derivative of the function. To perform this task, the automatic differentiation program  112  implements a variety of methods, which are described in greater detail below. 
     Method of Operating the Automatic Differentiation System 
       FIG. 2  shows a flow diagram for a method  200  for operating the automatic differentiation system. In the description of these method, statements that some task, calculation, or function is performed refers to a processor (e.g., the processor  102  of the automatic differentiation system  100 ) executing programmed instructions stored in non-transitory computer readable storage media (e.g., the memory  104  of the automatic differentiation system  100 ) operatively connected to the processor to manipulate data or to operate one or more components of the automatic differentiation system  100  to perform the task or function. Additionally, the steps of the methods may be performed in any feasible chronological order, regardless of the order shown in the figures or the order in which the steps are described. 
     The method  200  begins with a step of storing an automatic differentiation program, the automatic differentiation program including a program construct   that implements the mathematical derivative operator (block  210 ). Particularly, the memory  104  stores the automatic differentiation program  112 , which includes the derivative operator program construct  that implements mathematical derivative operator  . As used herein “program construct” refers to fragment of program code that can be invoked by a syntax that is identifiable by a compiler or interpreter. The automatic differentiation program  112  may, of course, include other program constructs that implement other operations. 
     The method  200  continues with a step of receiving first program code defining a function ƒ, the function ƒ being a higher-order function which takes an argument x as input (block  230 ). Particularly, the processor  102  is configured to receive first program code that defines a first mathematical function ƒ. As used herein, the term “function” refers to a concrete function that is implemented by computer program code. Accordingly, the first mathematical function ƒ and the first program code that implements or defines the first mathematical function ƒ can be considered essentially one and the same. 
     The first program code is configured to implement the first mathematical function ƒ to take an argument x as input and output a result based on the argument x. As used herein an “argument” or “independent variable” of a function refers to a value that is provided as an input to a function in order to obtain the function&#39;s result. In some cases, an argument may be a function which takes another argument as input to obtain the result. In at least one embodiment, the first program code is written in a functional programing language. 
     Additionally, the first mathematical function ƒ defined by the first program code is a higher-order function. As used herein a “higher-order” function refers a function that takes one or more functions as arguments and/or returns a function as its result. In other words, a higher-order has a range that is function and/or a domain that is a function. In contrast, as used herein a “first-order” function refers to a function that does not take a function as an argument and does not return a function as its result. In other words, a first-order function is a function having a non-function range and a non-function domain. Accordingly, the first mathematical function ƒ defined by the first program code is configured to take the argument x as input and output a function g which take an argument y. 
     In at least one embodiment, the memory  104  further stores some program P which has some arbitrary purpose and which requires the determination of a derivative of a mathematical function that is defined by a fragment of program code. The program P may utilize the automatic differentiation program  112  to determine the derivative of the mathematical function. In particular, the program P may invoke the derivative operator program construct  , providing the mathematical function that is defined by a fragment of program code as an argument. In the method  200 , the first mathematical function ƒ and the first program code is the fragment of program code that is provided as an argument to this top-level invocation the derivative operator program construct  . 
     In at least some embodiments, in addition to the first program code that defines a first mathematical function ƒ, the processor  102  further receives a value for the argument x. Particularly, in at least some cases, the program P may invoke the derivative operator program construct  , providing (i) the first program code that defines a first mathematical function ƒ and (ii) a value for the argument x, as arguments of the top-level invocation of the derivative operator program construct  . 
     With continued reference to  FIG. 2 , the method  200  continues with a step of determining, using the automatic differentiation program, a function ƒ′ which is the derivative of the function ƒ, the determination of the function ƒ′ including multiple executions of the program construct  , each execution being distinguished from one another using unique tags (block  250 ). Particularly, the processor  102  is configured to execute the automatic differentiation program  112  to determine a second mathematical function ƒ′, which is the derivative of the first mathematical function ƒ (i.e., ƒ′ x= ƒ x). The processor  102  determines the second mathematical function ƒ′ by applying the derivative operator program construct   of the automatic differentiation program  112  to the first mathematical function ƒ defined by the received first program code. 
     As used herein, “determining” the second mathematical function ƒ′ includes automatic differentiation using operator overloading techniques, source transformation techniques, or any other equivalent technique. Particularly, in some embodiments, the processor  102  utilizes source transformation techniques to generate second program code that defines the second mathematical function ƒ′, and which can be executed to evaluate the second mathematical function ƒ′ at particular values for the argument x. Additionally, in some embodiments, the processor  102  utilizes operator overloading techniques to evaluate the second mathematical function ƒ′ at particular values for the argument x using the first program code that defines the first mathematical function ƒ in conjunction with the derivative operator program construct  . 
     In the particular case that the first mathematical function ƒ is a higher-order function, the processor  102  must execute the derivative operator program construct   multiple times to determine the second mathematical function ƒ′. As used herein, an “execution” the derivative operator program construct  , refers to any distinct derivate calculation including a top-level invocation of the derivative operator program construct  , as well as any nested calculations using the derivative operator program construct  . 
     The derivative operator program construct   may implement any type of automatic differentiation. In some embodiments, the derivative operator program construct   implements forward mode AD. In some embodiments, the derivative operator program construct   implements reverse mode AD, which may utilize checkpointing. In some embodiments, the derivative operator program construct   implements a hybrid mode AD. 
     In some embodiments, the processor  102  advantageously utilizes unique tags to distinguish between each particular execution of the derivative operator program construct   or, in other words, each distinct derivative calculation. This can be accomplished by either of the first and second methods described above—that is, by eta expansion or tag substitution. 
     In some embodiments, the processor  102  determines the second mathematical function ƒ′ such that certain operations are performed at run-time and certain other operations are performed at compile-time. In this way, if feasible, any of the operations or steps discussed herein can be performed at run-time or at compile-time. Moreover, the processor  102  may determine whether a particular operation should be performed at run-time or at compile-time in a dynamic and intelligent manner. 
     In the determination of the second mathematical function ƒ′, the processor  102  implements three steps (which correspond essentially to the steps 1-3, described above) for each respective execution of the derivative operator program construct  . First, the processor  102  forms a dual number x+1ε 1  with the argument x as primal, with one as tangent, and with an infinitesimal number ε 1  having a unique tag (i.e., the subscript of the infinitesimal number ε 1 ) that uniquely associates it with the respective execution of the derivative operator program construct. As used herein, the term “tag” refers to any identifier, numeral, dependent data type, or other distinguishing feature that uniquely associates a particular infinitesimal number with a particular execution of the derivative operator program construct  . 
     Next, the processor  102  determines a result ƒ(x+1ε 1 ) by applying the first mathematical function ƒ with the respective dual number x+1ε 1  as input. Finally, the processor  102  extracts a tangent from the result ƒ(x+1ε 1 ) with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct, i.e. tg ε 1  ƒ(x+1ε 1 ). These steps are performed for each particular execution of the derivative operator program construct  . 
     As used herein a “dual number” refers to a number defined by three components: a primal, a tangent, and infinitesimal number ∈, where ∈ is an abstract number having the properties ∈ 2 =0 and ∈≠0. The abstract number e may be referred to herein as an “infinitesimal number” or a “nilpotent number.” In rectangular or Cartesian form, a dual number takes the form of a+b∈. However, it will be appreciated that mathematically equivalent forms may also be utilized, such as a polar coordinate form, a data triple, or data tuple having a dependent data type. The coefficient b of the infinitesimal number ∈ is referred to herein as the “tangent” or “perturbation” of the dual number. Similarly, the value a is referred to herein as the “primal” of the dual number. As used herein, “forming” a dual number with a primal a, a tangent b, and infinitesimal number ∈ refers to forming the dual number as a+b∈ or in any mathematically or computationally equivalent form. 
     It will be appreciated that a number may include multiple infinitesimal numbers having multiple different tags (e.g., ε 1 , ε 2 , . . . , ε n ) and can thus be considered a dual number with respect to each of the differently tagged infinitesimal numbers. For example, a number that includes both ε 1  and ε 2  may be manipulated into the form a+bε 1 , as well as into the form c+dε 2 . In this example, the value a is the primal with respect to ε 1  and the value b is the tangent with respect to ε 1 . Likewise, the value c is the primal with respect to ε 2  and the value d is the tangent with respect to ε 2 . Accordingly, as used herein a tangent “with respect to” a particular infinitesimal number ε (notated as tg ε, above) refers to the value b when the dual number is manipulated into the form a+bε. Likewise, as used herein a primal “with respect to” a particular infinitesimal number ε refers to the value a when the dual number is manipulated into the form a+bε. Thus, the primal a and tangent b with respect to a particular infinitesimal number (e.g., ε 1 ) may include another differently tagged infinitesimal number (e.g., ε 2 ), but will not include the particular infinitesimal number (e.g., ε 1 ) with respect to which the primal and tangent where extracted. 
     In embodiments utilizing eta expansion (first method), the processor  102  eta-expands the first mathematical function until all of the multiple executions of the derivative operator program construct   will result in non-function-containing results, in accordance with equations (18a) and (18b). This eta expansion is performed prior to the forming of the dual number x+1ε 1 , the determining of the result ƒ(x+1ε 1 ), and the extracting of the tangent tg ε 1  ƒ(x+1ε 1 ), for each of the multiple executions of the derivative operator program construct  . The processor  102  determines the second mathematical function ƒ′ as a result of the multiple executions of the derivative operator program construct   in the eta expansion of the first mathematical function ƒ according to equation. In some embodiments, prior to eta expansion, the processor  102  checks whether collision or confusion between distinct infinitesimal numbers is possible, and the eta expansion is terminated or omitted if collision or confusion is impossible. 
     In embodiments utilizing tag substitution (second method), the processor  102  determines the second mathematical function ƒ′ as being equal to the tangent tg ε 1  ƒ(x+1ε 1 ) in the outermost calculation using the derivative operator program construct  , in accordance with equation (8d), i.e.   ƒ x   fresh ε in tg ε ƒ(x+1ε). It will be appreciated that, in embodiments utilizing tag substitution, calculations using the derivative operator program construct   are nested and the outermost calculation refers to the final completed calculation using the using the derivative operator program construct  . 
     Regardless of whether eta expansion (first method) or tag substitution (second method) is utilized, for each respective execution of the derivative operator program construct  , the processor  102  extracts the tangent tg ε 1  ƒ(x+1ε 1 ) from the result ƒ(x+1ε 1 ) according to equations (8a), (8b), and (8c). Particularly, in response to the result ƒ(x+1ε 1 ) defining a real number, the processor  102  determines the tangent tg ε 1  ƒ(x+1ε 1 ) as zero in accordance with equation (8a), i.e. tg ε 1  a 0. Additionally, in response to the result ƒ(x+1ε 1 ) defining a dual number a+bε 1  having only the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct, the processor  102  determines the tangent tg ε 1  ƒ(x+1ε 1 ) as a tangent b of the dual number a+bε 1  with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct, in accordance with equation (8b), i.e. tg ε 1  (a+bε 1 ) b. Finally, in response to the result ƒ(x+1ε 1 ) defining a dual number c+dε 2  having an infinitesimal number ε 2  associated with a different execution of the derivative operator program construct, the processor  102  determines the first tangent tg ε 1  ƒ(x+1ε 1 ) according to equation (8c). Particularly, the processor  102  extracts a tangent tg ε 1  c with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct from a primal c of the dual number c+dε 2  with respect to the infinitesimal number ε 2  associated with the different execution of the derivative operator program construct. The processor  102  extracts extracting a tangent tg ε 1  d with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct from a tangent d of the dual number c+dε 2  with respect to the infinitesimal number ε 2  associated with the different execution of the derivative operator program construct. The processor  102  determines the first tangent tg ε 1  ƒ(x+1ε 1 ) as a dual number formed with the tangent tg ε 1  c as primal, with the tangent tg ε 1  d as tangent, and with the infinitesimal number ε 2  associated with the different execution of the derivative operator program construct, in accordance with equation (8c), e.g. tg ε 1  (c+dε 2 ) (tg ε 1  c)+(tg ε 1  d)ε 2 . 
     In the embodiments utilizing tag substitution (second method), for each respective execution of the derivative operator program construct  , the processor  102  may also extract the tangent tg ε 1  ƒ(x+1ε 1 ) from the result ƒ(x+1ε 1 ) according to equation (23). Particularly, in response to the result ƒ(x+1ε 1 ) defining a mathematical function  g  which takes an argument y as input, the processor  102  extracts the tangent tg ε 1  ƒ(x+1ε 1 ) according to equation (23). First, the processor  102  determines a result [ε/ε 1 ]y by substituting, in the argument y, the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct with a temporary infinitesimal number ε having a temporary unique tag. Next, the processor  102  determines a result  g ∘[ε/ε 1 ]y by applying the mathematical function g to the result [ε/ε 1 ]y. Next, the processor  102  extracts a tangent (tg ε 1 )∘ g ∘[ε/ε 1 ]y with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct from the result  g ∘[ε/ε 1 ]y. Finally, the processor  102  determines the tangent tg ε 1  ƒ(x+1ε 1 ) by substituting, in the tangent (tg ε 1 )∘ g ∘[ε/ε 1 ]y, the temporary infinitesimal number having the temporary unique tag with the infinitesimal number associated with the respective execution of the derivative operator program construct, according to equation (23), i.e. tg ε 1   g y≙fresh ε in ([ε 1 /ε]∘(tg ε 1 )∘ g ∘[ε/ε 1 ])y. 
     It will be appreciated that the determining the result  g ∘[ε/ε 1 ]y may involve a nested execution of the derivative operator program construct  . The nested execution of the derivative operator program construct includes forming a dual number with a further infinitesimal number ε 2  having a further unique tag. In this way, the tag substitution embodiments involved nested calculations using the derivative operator program construct  , in which each calculation is provided a fresh unique tag. 
     In the embodiments utilizing tag substitution (second method), for each performance of the substitution operation [ε 1 /ε 2 ], the processor  102  determines a respective substitution output by substituting, in a respective substitution input, infinitesimal numbers ε 2  having a second unique tag with infinitesimal numbers ε 1  having a first unique tag according to equations (24a), (24b), (24c), and (24d). Particularly, in response to the respective substitution input defining a real number, the processor  102  determines the respective substitution output as being equal to the respective substitution input, in accordance with equation (24a), i.e. [ε 1 /ε 2 ]a≙a. In response to the respective substitution input defining a dual number a+bε 2  having only infinitesimal numbers ε 2  having the second unique tag, the processor  102  determines the respective substitution output by substituting, in the dual number a+bε 2 , infinitesimal numbers ε 2  having the second unique tag with infinitesimal numbers ε 1  having the first unique tag, in accordance with equation (24b), i.e. [ε 1 /ε 2 ](a+bε 2 ) a+bε 1 . 
     In response to the respective substitution input defining a dual number a+bε having infinitesimal numbers e having a further unique tag, the processor  102  determines the respective substitution output according to equation (24c). Particularly, the processor  102  determines a result [ε 1 /ε 2 ]a by substituting infinitesimal numbers ε 2  having the second unique tag with infinitesimal numbers ε 1  having the first unique tag, in a primal a of the dual number a+bε with respect to the infinitesimal numbers e having the further unique tag. Next, the processor  102  determines a result [ε 1 /ε 2 ]b by substituting infinitesimal numbers ε 2  having the second unique tag with infinitesimal numbers ε 1  having the first unique tag, in a tangent b of the dual number a+bε with respect to the infinitesimal numbers e having the further unique tag. Finally, the processor  102  determines the respective substitution output as a dual number formed with the result [ε 1 /ε 2 ]a as primal, with the result [ε 1 /ε 2 ]b as tangent, and with the infinitesimal number e having the further unique tag, in accordance with equation (24c), i.e. [ε 1 /ε 2 ](a+bε) ([ε 1 /ε 2 ]a)+([ε 1 /ε 2 ]b)ε. 
     In response to the respective substitution input defining a mathematical function g which takes an argument y as input, the processor  102  determines the respective substitution output according to equation (24d). Particularly, the processor  102  determines a result [ε/ε 2 ]y by substituting infinitesimal numbers ε 2  having the second unique tag with infinitesimal numbers e having a temporary unique tag in the argument y. Next, the processor  102  determines a result  g ∘[ε/ε 2 ]y as derivative of the result [ε/ε 2 ]y. Next, the processor  102  determines a result [ε 1 /ε 2 ]∘ g ∘[ε/ε 2 ]y by substituting infinitesimal numbers ε 2  having the second unique tag with infinitesimal numbers ε 1  having the first unique tag in the result  g ∘[ε/ε 2 ]y. Finally, the processor  102  determines the respective substitution output by substituting infinitesimal numbers ε having the temporary unique tag with infinitesimal numbers having the second unique tag ε 2  in the result, in accordance with equation (23), i.e. [ε 1 /ε 2 ] g y fresh ε in ([ε 2 /ε]∘[ε 1 /ε 2 ]∘ g ∘[ε/ε 2 ])y. 
     It should be appreciated that, depending on the programming languages used, the substitution operation described above can equivalently be performed using other mechanisms by wrapping the derivative calculations in prologue and epilogue code that creates new distinguishing features and substitutes them for the distinguishing feature of the wrapped procedure in its argument before invoking the wrapped procedure, and reverses the process for the output of the wrapped procedure. 
     In some embodiments, the first program code may define a first mathematical function ƒ whose range and domain are both functions or, in other words, the first mathematical ƒ takes arguments x and x′ as inputs, which may be themselves functions of an argument y. In these cases, the processor  102  is configured to utilize a derivative operator program construct   instead of the derivative operator program construct  . 
     In the determination of the second mathematical function ƒ′, the processor  102  implements three steps for each respective execution of the derivative operator program construct  . First, the processor  102  forms a bundle bun ε 1  x x′ including with the arguments x and x′ with respect to an infinitesimal number ε 1  having a unique tag (i.e., the subscript of the infinitesimal number ε 1 ) that uniquely associates it with the respective execution of the derivative operator program construct. Next, the processor  102  determines a result ƒ(bun ε 1  x x′) by applying the first mathematical function ƒ with the bundle bun ε 1  x x′ as input. Finally, the processor  102  extracts a tangent from the result ƒ(bun ε 1  x x′) with respect to the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct, i.e. tg ε 1  ƒ(bun ε 1  x x′), in accordance with equations (32b) or (33). These steps are performed for each particular execution of the derivative operator program construct  . 
     As used herein, a “bundle” formed with a first argument and a second argument with respect to a particular infinitesimal number refers to a dual number formed with first argument as primal, with the second argument as tangent, and with the particular infinitesimal number. In the case that the first and second arguments are functions of a third argument, then the functions are applied to the third argument before forming the dual number. 
     In embodiments utilizing eta expansion (first method), the processor  102  eta-expands the first mathematical function until all of the multiple executions of the derivative operator program construct   will result in non-function-containing results, in accordance with equations (32a) and (32b). This eta expansion is performed prior to the forming of the bundle bun ε 1  x x′, the determining of the result ƒ(bun ε 1  x x′), and the extracting of the tangent tg ε 1  ƒ(bun ε 1  x x′), for each of the multiple executions of the derivative operator program construct  . The processor  102  determines the second mathematical function ƒ′ as a result of the multiple executions of the derivative operator program construct   in the eta expansion of the first mathematical function ƒ according to equation. 
     In embodiments utilizing tag substitution (second method), the processor  102  determines the second mathematical function ƒ′ as being equal to the tangent tg ε 1  ƒ(bun ε 1  x x′) in the outermost calculation using the derivative operator program construct  , in accordance with equation (33), i.e.   ƒ x x′ fresh ε in tg ε ƒ(bun ε 1  x x′). It will be appreciated that, in embodiments utilizing tag substitution, calculations using the derivative operator program construct   are nested and the outermost calculation refers to the final completed calculation using the using the derivative operator program construct  . 
     Regardless of whether eta expansion (first method) or tag substitution (second method) is utilized, for each respective execution of the derivative operator program construct  , in response to the arguments x and x′ not defining or containing functions, the processor  102  determines the bundle bun ε 1  x x′ as a dual number x+x′ε 1  formed with the first argument x as primal, with the second argument x′ as tangent, and with the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct  . However, as noted below, in at least some cases the arguments x and x′ are instead some functions of an argument y, i.e. x=(ƒ y) and x′=(ƒ′ y). 
     In the embodiments utilizing eta expansion (first method), for each respective execution of the derivative operator program construct  , the processor  102  also determines the bundle bun ε 1  x x′ according to equation (30b), in response to the arguments x and x′ defining some functions of a further argument y, i.e. x=(ƒ y) and x′=(ƒ′ y). Particularly, the processor  102  determines a result (ƒ y) by applying the mathematical function x=f with the argument y as input. Next, the processor  102  determines a result (ƒ′ y) by applying the mathematical function x′=ƒ′ with the argument y as input. Finally, the processor and (iii) determining the bundle bun ε 1  x x′=bun ε 1  ƒ ƒ′y as a dual number (ƒ′ y)+(ƒ′ y)ε 1  formed with the result x=(ƒ y) as primal, with the result x′=(ƒ′ y) as tangent, and with the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct  , i.e. bun ε 1  ƒ ƒ′y bun ε 1  (ƒy)(ƒ′y). 
     In the embodiments utilizing tag substitution (second method), for each respective execution of the derivative operator program construct  , the processor  102  may also determine the bundle bun ε 1  x x′ according to equation (31), in response to the arguments x and x′ defining some functions of a further argument y, i.e. x=(ƒ y) and x′=(ƒ′ y), Particularly, the processor  102  determines a result [ε/ε 1 ]y by substituting infinitesimal numbers ε 1  having the second unique tag with infinitesimal numbers ε having a temporary unique tag in the argument y. Next, the processor  102  determines a result (ƒ ([ε/ε 1 ]y)) by applying the mathematical function x=ƒ with the result [ε/ε 1 ]y as input. Next, the processor  102  determines a result (ƒ′([ε/ε 1 ]y)) by applying the mathematical function x′=ƒ′ with the result [ε/ε 1 ]y as input. Next, the processor  102  determines a bundle bun ε 1  (ƒ[ε/ε 1 ]y)(ƒ′[ε/ε 1 ]y)) as a dual number (ƒ ([ε/ε 1 ]y))+(ƒ′([ε/ε 1 ]y))ε 1  formed with the result (ƒ ([ε/ε 1 ]y)) as primal, with the result (ƒ′([ε/ε 1 ]y)) as tangent, and with the infinitesimal number ε 1  associated with the respective execution of the derivative operator program construct  . Finally, the processor  102  determines the bundle bun ε 1  x x′=bun ε 1  ƒ ƒ′y by substituting, in the bundle bun ε 1  (ƒ[ε/ε 1 ]y)(ƒ′[ε/ε 1 ]y)), the infinitesimal numbers e having a temporary unique tag with the infinitesimal numbers ε 1  associated with the respective execution of the derivative operator program construct  , i.e. bun ε 1  ƒ ƒ′y fresh ε in [ε/ε 1 ](bun ε 1  (ƒ[ε/ε 1 ]y)(ƒ′[ε/ε 1 ]y)). 
     With continued reference to  FIG. 2 , the method  200  continues with a step of returning, outputting, or storing at least one of (i) second program code defining function ƒ′, and (ii) a result of the function ƒ′ at one or more values for the argument x. (block  270 ). Particularly, the processor  102  is configured to return, output, or store at least one of (i) second program code defining function ƒ′, and (ii) a result of the function ƒ′ at one or more values for the argument x. As noted above, in some embodiments a program P, which has some arbitrary purpose, may invoke the derivative operator program construct  , providing the first program code defining the first mathematical function ƒ as an argument of the derivative operator program construct  . In some cases, the program P provides one or more particular values for the argument x as a further argument of the derivative operator program construct  . The processor  102  returns to the program P at least one of (i) second program code defining second mathematical function ƒ′, and (ii) a result of the second mathematical function ƒ′ at the one or more values for the argument x. In practice, this may include storing the returns the memory  104  or otherwise outputting the returns such that they are provided to and accessibly by the program P. 
     In embodiments utilizing source transformation, the processor  102  may generate and return the second program code defining second mathematical function ƒ′, which can then be executed in the program P to evaluate the second mathematical function ƒ′ at one or more values for the argument x. However, in embodiments utilizing operator overloading, the processor  102  may simply evaluate and return results of the second mathematical function ƒ′ at the one or more values for the argument x. 
     Exemplary Implementations of the Automatic Differentiation Program 
       FIGS. 3A-3B  and  FIGS. 4A-4C  show program code for exemplary minimal implementations of the automatic differentiation program  112 . These exemplary minimal implementations are not intended as a full practical implementation but rather has the expository purpose of explaining the ideas presented in this disclosure. 
       FIGS. 3A-3B  show program code  300  for exemplary minimal implementations of the automatic differentiation program  112  which implements the first method described above that utilizes the eta expansion technique. The program code  300  includes a program code fragment  302 , which implements the tg operation using equations (8a), (8b), and (8c). The program code  300  further includes a program code fragment  304 , which implements the fresh E operation. The program code  300  further includes a program code fragment  306 , which implements the   operation that materializes D using equations (18a) and (18b). The program code  300  further includes a program code fragment  308 , which implements the bun operation using equations (30a) and (30b). The program code  300  further includes a program code fragment  310 , which implements the   operation that materializes   using equations (32a) and (32b). 
       FIGS. 4A-4C  show program code  400  for exemplary minimal implementations of the automatic differentiation program  112  which implements the second method described above that utilizes the tag substitution technique. The program code  400  includes a program code fragment  402 , which implements the tg operation using equations (8a), (8b), (8c), and (23). The program code  400  further includes a program code fragment  404 , which implements the [ε 1 /ε 2 ] operation using equations (24a), (24b), (24c), and (24d). The program code  400  further includes a program code fragment  406 , which implements the fresh E operation. The program code  400  further includes a program code fragment  408 , which implements the   operation that materializes   using equations (8d). The program code  400  further includes a program code fragment  410 , which implements the bun operation using equations (30a) and (31). The program code  400  further includes a program code fragment  412 , which implements the   operation that materializes   using equation (33). 
     While the disclosure has been illustrated and described in detail in the drawings and foregoing description, the same should be considered as illustrative and not restrictive in character. It is understood that only the preferred embodiments have been presented and that all changes, modifications and further applications that come within the spirit of the disclosure are desired to be protected.