\lmcsdoi

18141 \lmcsheadingLABEL:LastPageJan. 19, 2021Mar. 22, 2022

Higher Order Automatic Differentiation of Higher Order Functions

Mathieu Huot\rsupera , Sam Staton\rsupera and Matthijs Vákár\rsuperb University of Oxford Utrecht University

Abstract.

We present semantic correctness proofs of automatic differentiation (AD). We consider a forward-mode AD method on a higher order language with algebraic data types, and we characterise it as the unique structure preserving macro given a choice of derivatives for basic operations. We describe a rich semantics for differentiable programming, based on diffeological spaces. We show that it interprets our language, and we phrase what it means for the AD method to be correct with respect to this semantics. We show that our characterisation of AD gives rise to an elegant semantic proof of its correctness based on a gluing construction on diffeological spaces. We explain how this is, in essence, a logical relations argument. Throughout, we show how the analysis extends to AD methods for computing higher order derivatives using a Taylor approximation.

Key words and phrases:

automatic differentiation, software correctness, denotational semantics

The authors contributed equally to this work.

1. Introduction

Automatic differentiation (AD), loosely speaking, is the process of taking a program describing a function, and constructing the derivative of that function by applying the chain rule across the program code. As gradients play a central role in many aspects of machine learning, so too do automatic differentiation systems such as TensorFlow [AAB⁺16], PyTorch [PGC⁺17] or Stan [CHB⁺15].

Figure 1. Overview of semantics/correctness of AD.

Differentiation has a well-developed mathematical theory in terms of differential geometry. The aim of this paper is to formalize this connection between differential geometry and the syntactic operations of AD, particularly for AD methods that calculate higher order derivatives. In this way we achieve two things: (1) a compositional, denotational understanding of differentiable programming and AD; (2) an explanation of the correctness of AD.

This intuitive correspondence (summarized in Fig. 1) is in fact rather complicated. In this paper, we focus on resolving the following problem: higher order functions play a key role in programming, and yet they have no counterpart in traditional differential geometry. Moreover, we resolve this problem while retaining the compositionality of denotational semantics.

1.0.1. Higher order functions and differentiation.

A major application of higher order functions is to support disciplined code reuse. Code reuse is particularly acute in machine learning. For example, a multi-layer neural network might be built of millions of near-identical neurons, as follows.

\begin{array}[]{ll}\begin{aligned} &\mathrm{neuron}_{n}:\boldsymbol{(}\mathbf{% real}^{n}\boldsymbol{\mathop{*}}\boldsymbol{(}\mathbf{real}^{n}\boldsymbol{% \mathop{*}}\mathbf{real}\boldsymbol{)}\boldsymbol{)}\to\mathbf{real}\\ &\mathrm{neuron}_{n}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda\langle x% ,\langle w,b\rangle\rangle.\,\varsigma(w\cdot x+b)\\ &\mathrm{layer}_{n}:(\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}P% \boldsymbol{)}\to{\tau}_{2})\to\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}% P^{n}\boldsymbol{)}\to{\tau}_{2}^{n}\\ &\mathrm{layer}_{n}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda f.\,% \lambda\langle x,\langle p_{1},\dots,p_{n}\rangle\rangle.\,\langle f\langle x,% p_{1}\rangle,\dots,f\langle x,p_{n}\rangle\rangle\\ &\mathrm{comp}:\boldsymbol{(}(\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}P% \boldsymbol{)}\to{\tau}_{2})\boldsymbol{\mathop{*}}(\boldsymbol{(}{\tau}_{2}% \boldsymbol{\mathop{*}}Q\boldsymbol{)}\to{\tau}_{3})\boldsymbol{)}\to% \boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}\boldsymbol{(}P\boldsymbol{% \mathop{*}}Q\boldsymbol{)}\boldsymbol{)}\to{\tau}_{3}\\ &\mathrm{comp}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda\langle f,g% \rangle.\,\lambda\langle x,(p,q)\rangle.\,g\langle f\langle x,p\rangle,q% \rangle\end{aligned}&\raisebox{-22.76219pt}[0.0pt]{ \leavevmode\hbox to106.35% pt{\vbox to89.59pt{\pgfpicture\makeatletter\hbox{\hskip 34.9209pt\lower-24.462% 46pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}% \pgfsys@invoke{ }\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\hbox to0.0pt{% \hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to% 0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{% \hbox to0.0pt{\hbox to0.0pt{\hbox to0.0pt{ {}{}{} \hss} {}{}{} \hss} {}{}{} { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{}{0.0pt}\pgfsys@invoke{ }% \definecolor[named]{tikz@color}{rgb}{.5,.5,.5}\definecolor[named]{.}{rgb}{% .5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}% \pgfsys@color@gray@stroke{.5}\pgfsys@invoke{ }\pgfsys@color@gray@fill{.5}% \pgfsys@invoke{ }\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{ }% \pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\color[rgb]{0.75,0.75,0.75}% \definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}{}% \pgfsys@moveto{0.0pt}{-0.00737pt}\pgfsys@lineto{71.13188pt}{-0.00737pt}% \pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{}{0.0pt}\pgfsys@invoke{ }% \definecolor[named]{tikz@color}{rgb}{.5,.5,.5}\definecolor[named]{.}{rgb}{% .5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}% \pgfsys@color@gray@stroke{.5}\pgfsys@invoke{ }\pgfsys@color@gray@fill{.5}% \pgfsys@invoke{ }\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{ }% \pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\color[rgb]{0.75,0.75,0.75}% \definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}{}% \pgfsys@moveto{0.0pt}{30.77252pt}\pgfsys@lineto{71.13188pt}{30.77252pt}% \pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setdash{}{0.0pt}\pgfsys@invoke{ }% \definecolor[named]{tikz@color}{rgb}{.5,.5,.5}\definecolor[named]{.}{rgb}{% .5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}% \pgfsys@color@gray@stroke{.5}\pgfsys@invoke{ }\pgfsys@color@gray@fill{.5}% \pgfsys@invoke{ }\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{ }% \pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\color[rgb]{0.75,0.75,0.75}% \definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}{}% \pgfsys@moveto{0.0pt}{61.55244pt}\pgfsys@lineto{71.13188pt}{61.55244pt}% \pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{15.9682pt}{-4.89998pt}\pgfsys@lineto{15.9682% pt}{-3.49998pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{35.92845pt}{-4.89998pt}\pgfsys@lineto{35.928% 45pt}{-3.49998pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{55.8887pt}{-4.89998pt}\pgfsys@lineto{55.8887% pt}{-3.49998pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{-0.00737pt}\pgfsys@lineto{-2.5pt}{-0% .00737pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{30.77252pt}\pgfsys@lineto{-2.5pt}{30% .77252pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {}{{}}{} {}{}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{61.55244pt}\pgfsys@lineto{-2.5pt}{61% .55244pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }{}\pgfsys@moveto{0.% 0pt}{-3.49998pt}\pgfsys@lineto{71.13188pt}{-3.49998pt}\pgfsys@stroke% \pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\color[rgb]{% 0.75,0.75,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}% \pgfsys@rectcap\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{7% 1.13188pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\color[rgb]{% 0.75,0.75,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}% \pgfsys@rectcap\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{61.54305pt}% \pgfsys@lineto{71.13188pt}{61.54305pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }{}\pgfsys@moveto{-2% .5pt}{0.0pt}\pgfsys@lineto{-2.5pt}{61.54305pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\color[rgb]{% 0.75,0.75,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}% \pgfsys@rectcap\pgfsys@invoke{ }{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{0% .0pt}{61.54305pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } \pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\color[rgb]{% 0.75,0.75,0.75}\definecolor[named]{pgfstrokecolor}{rgb}{0.75,0.75,0.75}% \pgfsys@color@gray@stroke{0.75}\pgfsys@invoke{ }\pgfsys@color@gray@fill{0.75}% \pgfsys@invoke{ }\definecolor{pgffillcolor}{rgb}{0.75,0.75,0.75}% \pgfsys@rectcap\pgfsys@invoke{ }{}\pgfsys@moveto{71.13188pt}{0.0pt}% \pgfsys@lineto{71.13188pt}{61.54305pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \hss}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hbox to0.0pt{% \pgfsys@beginscope\pgfsys@invoke{ }{{}}{{}} {} {}{}{}{}{}{}{}{}{}{} {} {} {}{} {} {} {}{}{}{} {} {} {} {}{}{}{} {}{}{}{} {} {} {} {}{}{}{} {} {} {} {}{} {}{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} { {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} } { } { } {} {} {} { {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} { } {} } \pgfsys@beginscope\pgfsys@invoke{ } \pgfsys@beginscope\pgfsys@invoke{ }{} \pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\definecolor[named]{tikz@color}{% rgb}{0,0,0}\definecolor[named]{.}{rgb}{0,0,0}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ } {}{}{}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}}{}{} {}{{}{}{}} {}{}{}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@curveto{0.0pt}{0.0pt}{0.52441pt}{0% .00133pt}{0.72583pt}{0.00185pt}\pgfsys@curveto{0.92725pt}{0.00235pt}{1.25024pt% }{0.00319pt}{1.45166pt}{0.0037pt}\pgfsys@curveto{1.65308pt}{0.0042pt}{1.97609% pt}{0.00485pt}{2.1775pt}{0.00552pt}\pgfsys@curveto{2.37892pt}{0.0062pt}{2.7019% 2pt}{0.00775pt}{2.90334pt}{0.0086pt}\pgfsys@curveto{3.10475pt}{0.00946pt}{3.42% 776pt}{0.01076pt}{3.62917pt}{0.01169pt}\pgfsys@curveto{3.83058pt}{0.01262pt}{4% .15358pt}{0.01428pt}{4.355pt}{0.01538pt}\pgfsys@curveto{4.55641pt}{0.01648pt}{% 4.87943pt}{0.01833pt}{5.08084pt}{0.01968pt}\pgfsys@curveto{5.28226pt}{0.02104% pt}{5.60526pt}{0.02362pt}{5.80667pt}{0.02524pt}\pgfsys@curveto{6.00809pt}{0.02% 684pt}{6.33109pt}{0.02943pt}{6.5325pt}{0.03139pt}\pgfsys@curveto{6.73392pt}{0.% 03334pt}{7.05693pt}{0.037pt}{7.25835pt}{0.03938pt}\pgfsys@curveto{7.45976pt}{0% .04176pt}{7.78276pt}{0.0458pt}{7.98418pt}{0.04863pt}\pgfsys@curveto{8.1856pt}{% 0.05144pt}{8.50859pt}{0.05629pt}{8.71pt}{0.0597pt}\pgfsys@curveto{8.91142pt}{0% .06311pt}{9.23444pt}{0.06915pt}{9.43585pt}{0.07324pt}\pgfsys@curveto{9.63727pt% }{0.07733pt}{9.96027pt}{0.0844pt}{10.16168pt}{0.08925pt}\pgfsys@curveto{10.363% 1pt}{0.09412pt}{10.6861pt}{0.10246pt}{10.88751pt}{0.10834pt}\pgfsys@curveto{11% .08893pt}{0.11423pt}{11.41193pt}{0.12465pt}{11.61334pt}{0.13173pt}% \pgfsys@curveto{11.81476pt}{0.13881pt}{12.13777pt}{0.15099pt}{12.33919pt}{0.15% 942pt}\pgfsys@curveto{12.5406pt}{0.16788pt}{12.8636pt}{0.1826pt}{13.06502pt}{0% .19267pt}\pgfsys@curveto{13.26643pt}{0.20274pt}{13.58943pt}{0.22003pt}{13.7908% 5pt}{0.23207pt}\pgfsys@curveto{13.99226pt}{0.24411pt}{14.31528pt}{0.26495pt}{1% 4.5167pt}{0.27948pt}\pgfsys@curveto{14.71811pt}{0.29399pt}{15.0411pt}{0.3194pt% }{15.24252pt}{0.33673pt}\pgfsys@curveto{15.44394pt}{0.35406pt}{15.76694pt}{0.3% 8368pt}{15.96835pt}{0.40445pt}\pgfsys@curveto{16.16977pt}{0.4252pt}{16.49278pt% }{0.46147pt}{16.6942pt}{0.48631pt}\pgfsys@curveto{16.89561pt}{0.51115pt}{17.21% 861pt}{0.55394pt}{17.42003pt}{0.58357pt}\pgfsys@curveto{17.62144pt}{0.6132pt}{% 17.94444pt}{0.66449pt}{18.14586pt}{0.69992pt}\pgfsys@curveto{18.34727pt}{0.735% 37pt}{18.67029pt}{0.79669pt}{18.87169pt}{0.83905pt}\pgfsys@curveto{19.0731pt}{% 0.88141pt}{19.3961pt}{0.9547pt}{19.59752pt}{1.00526pt}\pgfsys@curveto{19.79893% pt}{1.05583pt}{20.12195pt}{1.14336pt}{20.32336pt}{1.20349pt}\pgfsys@curveto{20% .52478pt}{1.26361pt}{20.84778pt}{1.36708pt}{21.0492pt}{1.43864pt}% \pgfsys@curveto{21.25061pt}{1.51022pt}{21.57361pt}{1.6343pt}{21.77502pt}{1.719% 36pt}\pgfsys@curveto{21.97644pt}{1.80443pt}{22.29945pt}{1.95091pt}{22.50087pt}% {2.05177pt}\pgfsys@curveto{22.70229pt}{2.15263pt}{23.02528pt}{2.32697pt}{23.22% 67pt}{2.44637pt}\pgfsys@curveto{23.42812pt}{2.56578pt}{23.75111pt}{2.77155pt}{% 23.95253pt}{2.91238pt}\pgfsys@curveto{24.15395pt}{3.05322pt}{24.47696pt}{3.295% 72pt}{24.67838pt}{3.46149pt}\pgfsys@curveto{24.87979pt}{3.62727pt}{25.20279pt}% {3.91287pt}{25.4042pt}{4.10727pt}\pgfsys@curveto{25.60562pt}{4.30165pt}{25.928% 62pt}{4.63567pt}{26.13004pt}{4.86261pt}\pgfsys@curveto{26.33145pt}{5.08955pt}{% 26.65446pt}{5.47916pt}{26.85588pt}{5.7429pt}\pgfsys@curveto{27.0573pt}{6.00665% pt}{27.3803pt}{6.45882pt}{27.58171pt}{6.76357pt}\pgfsys@curveto{27.78313pt}{7.% 06831pt}{28.10612pt}{7.58952pt}{28.30754pt}{7.93936pt}\pgfsys@curveto{28.50896% pt}{8.28922pt}{28.83195pt}{8.88652pt}{29.03337pt}{9.28505pt}\pgfsys@curveto{29% .23479pt}{9.6836pt}{29.5578pt}{10.36154pt}{29.75922pt}{10.81175pt}% \pgfsys@curveto{29.96063pt}{11.26196pt}{30.28363pt}{12.02603pt}{30.48505pt}{12% .52988pt}\pgfsys@curveto{30.68646pt}{13.03372pt}{31.00946pt}{13.88542pt}{31.21% 088pt}{14.44316pt}\pgfsys@curveto{31.4123pt}{15.00092pt}{31.7353pt}{15.93964pt% }{31.93672pt}{16.54973pt}\pgfsys@curveto{32.13814pt}{17.15984pt}{32.46114pt}{1% 8.1816pt}{32.66255pt}{18.84038pt}\pgfsys@curveto{32.86397pt}{19.49918pt}{33.18% 697pt}{20.59636pt}{33.38838pt}{21.29785pt}\pgfsys@curveto{33.5898pt}{21.99934% pt}{33.91281pt}{23.16028pt}{34.11423pt}{23.89629pt}\pgfsys@curveto{34.31564pt}% {24.6323pt}{34.63864pt}{25.84221pt}{34.84006pt}{26.60246pt}\pgfsys@curveto{35.% 04146pt}{27.36272pt}{35.36447pt}{28.603pt}{35.56589pt}{29.37573pt}% \pgfsys@curveto{35.7673pt}{30.14845pt}{36.0903pt}{31.39896pt}{36.29172pt}{32.1% 7177pt}\pgfsys@curveto{36.49313pt}{32.94458pt}{36.81615pt}{34.18541pt}{37.0175% 6pt}{34.94566pt}\pgfsys@curveto{37.21898pt}{35.70592pt}{37.54198pt}{36.9153pt}% {37.7434pt}{37.65121pt}\pgfsys@curveto{37.94481pt}{38.38715pt}{38.2678pt}{39.5% 4817pt}{38.46922pt}{40.24966pt}\pgfsys@curveto{38.67064pt}{40.95114pt}{38.9936% 5pt}{42.04834pt}{39.19507pt}{42.70714pt}\pgfsys@curveto{39.39648pt}{43.36592pt% }{39.71948pt}{44.38776pt}{39.9209pt}{44.99777pt}\pgfsys@curveto{40.12231pt}{45% .60779pt}{40.44531pt}{46.54607pt}{40.64673pt}{47.10373pt}\pgfsys@curveto{40.84% 814pt}{47.66139pt}{41.17114pt}{48.51324pt}{41.37256pt}{49.01701pt}% \pgfsys@curveto{41.57397pt}{49.52077pt}{41.89699pt}{50.28432pt}{42.0984pt}{50.% 73453pt}\pgfsys@curveto{42.29982pt}{51.18474pt}{42.62282pt}{51.8633pt}{42.8242% 3pt}{52.26183pt}\pgfsys@curveto{43.02565pt}{52.66037pt}{43.34866pt}{53.25714pt% }{43.55006pt}{53.6069pt}\pgfsys@curveto{43.75148pt}{53.95667pt}{44.07448pt}{54% .47803pt}{44.2759pt}{54.7827pt}\pgfsys@curveto{44.47731pt}{55.08737pt}{44.8003% 2pt}{55.53908pt}{45.00174pt}{55.80275pt}\pgfsys@curveto{45.20316pt}{56.0664pt}% {45.52615pt}{56.45612pt}{45.72757pt}{56.68306pt}\pgfsys@curveto{45.92899pt}{56% .91pt}{46.25198pt}{57.244pt}{46.4534pt}{57.4384pt}\pgfsys@curveto{46.65482pt}{% 57.6328pt}{46.97783pt}{57.91838pt}{47.17924pt}{58.08415pt}\pgfsys@curveto{47.3% 8066pt}{58.24994pt}{47.70366pt}{58.49243pt}{47.90508pt}{58.63327pt}% \pgfsys@curveto{48.10649pt}{58.77411pt}{48.42949pt}{58.97997pt}{48.6309pt}{59.% 09927pt}\pgfsys@curveto{48.83232pt}{59.2186pt}{49.15533pt}{59.3924pt}{49.35675% pt}{59.49326pt}\pgfsys@curveto{49.55817pt}{59.59413pt}{49.88116pt}{59.74124pt}% {50.08258pt}{59.82631pt}\pgfsys@curveto{50.284pt}{59.91138pt}{50.60701pt}{60.0% 3482pt}{50.80841pt}{60.1064pt}\pgfsys@curveto{51.00983pt}{60.17798pt}{51.33282% pt}{60.28204pt}{51.53424pt}{60.34218pt}\pgfsys@curveto{51.73566pt}{60.4023pt}{% 52.05867pt}{60.48929pt}{52.26009pt}{60.53978pt}\pgfsys@curveto{52.4615pt}{60.5% 9026pt}{52.7845pt}{60.66364pt}{52.98592pt}{60.706pt}\pgfsys@curveto{53.18733pt% }{60.74835pt}{53.51033pt}{60.80968pt}{53.71175pt}{60.84511pt}\pgfsys@curveto{5% 3.91316pt}{60.88055pt}{54.23618pt}{60.93184pt}{54.43759pt}{60.96146pt}% \pgfsys@curveto{54.639pt}{60.99109pt}{54.962pt}{61.03387pt}{55.16342pt}{61.058% 73pt}\pgfsys@curveto{55.36484pt}{61.08359pt}{55.68784pt}{61.11986pt}{55.88925% pt}{61.14061pt}\pgfsys@curveto{56.09067pt}{61.16136pt}{56.41368pt}{61.191pt}{5% 6.6151pt}{61.20833pt}\pgfsys@curveto{56.81651pt}{61.22566pt}{57.13951pt}{61.25% 105pt}{57.34093pt}{61.26556pt}\pgfsys@curveto{57.54234pt}{61.28008pt}{57.86534% pt}{61.30084pt}{58.06676pt}{61.31297pt}\pgfsys@curveto{58.26817pt}{61.32509pt}% {58.59117pt}{61.34291pt}{58.79259pt}{61.35298pt}\pgfsys@curveto{58.994pt}{61.3% 6305pt}{59.31702pt}{61.37724pt}{59.51843pt}{61.3856pt}\pgfsys@curveto{59.71985% pt}{61.39397pt}{60.04285pt}{61.40623pt}{60.24426pt}{61.41331pt}\pgfsys@curveto% {60.44568pt}{61.4204pt}{60.76868pt}{61.43082pt}{60.9701pt}{61.43669pt}% \pgfsys@curveto{61.17151pt}{61.44258pt}{61.49452pt}{61.45093pt}{61.69594pt}{61% .45578pt}\pgfsys@curveto{61.89735pt}{61.46065pt}{62.22035pt}{61.46771pt}{62.42% 177pt}{61.4718pt}\pgfsys@curveto{62.62318pt}{61.47589pt}{62.94618pt}{61.48193% pt}{63.1476pt}{61.48534pt}\pgfsys@curveto{63.34901pt}{61.48875pt}{63.67203pt}{% 61.4936pt}{63.87344pt}{61.49641pt}\pgfsys@curveto{64.07486pt}{61.49922pt}{64.3% 9786pt}{61.50328pt}{64.59927pt}{61.50566pt}\pgfsys@curveto{64.80069pt}{61.5080% 4pt}{65.12369pt}{61.5117pt}{65.3251pt}{61.51366pt}\pgfsys@curveto{65.52652pt}{% 61.51561pt}{65.84953pt}{61.51819pt}{66.05093pt}{61.5198pt}\pgfsys@curveto{66.2% 5235pt}{61.52142pt}{66.57535pt}{61.52399pt}{66.77676pt}{61.52534pt}% \pgfsys@curveto{66.97818pt}{61.52672pt}{67.3012pt}{61.52856pt}{67.50261pt}{61.% 52966pt}\pgfsys@curveto{67.70403pt}{61.53076pt}{68.02702pt}{61.53242pt}{68.228% 44pt}{61.53336pt}\pgfsys@curveto{68.42986pt}{61.53429pt}{68.75285pt}{61.53557% pt}{68.95427pt}{61.53642pt}\pgfsys@curveto{69.15569pt}{61.53728pt}{69.4787pt}{% 61.53883pt}{69.68011pt}{61.53952pt}\pgfsys@curveto{69.88153pt}{61.54019pt}{70.% 20453pt}{61.54085pt}{70.40594pt}{61.54137pt}\pgfsys@curveto{70.60736pt}{61.541% 87pt}{71.13177pt}{61.54321pt}{71.13177pt}{61.54321pt}\pgfsys@stroke% \pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@beginscope\pgfsys@invoke{ }{} \pgfsys@setlinewidth{0.6pt}\pgfsys@invoke{ }\definecolor[named]{tikz@color}{% rgb}{0,0,0}\definecolor[named]{.}{rgb}{0,0,0}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ } \pgfsys@setdash{}{0.0pt}\pgfsys@invoke{ }\definecolor[named]{tikz@color}{rgb}{% 0,0,0}\definecolor[named]{.}{rgb}{0,0,0}\definecolor[named]{pgfstrokecolor}{% rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ } {}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope { }{ {}{} } {\pgfsys@beginscope\pgfsys@invoke{ }{}{}\pgfsys@setlinewidth{0.4pt}% \pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}\pgfsys@setlinewidth{0.4pt}% \pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}\pgfsys@setlinewidth{0.4pt}% \pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} { {}{} }{ {}{} } {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{15.9682pt}{-4.89998pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{11.52377pt}{-12.72153pt}\pgfsys@invoke{ }\hbox{{\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$-5$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{35.92845pt}{-4.89998pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{33.92845pt}{-12.72153pt}\pgfsys@invoke{ }\hbox{{\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{55.8887pt}{-4.89998pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{53.8887pt}{-12.72153pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor% }{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$5$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {\pgfsys@beginscope\pgfsys@invoke{ }{}{}{}{}{}{}{}{} {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{-0.00737pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{{}{}}}% {{}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{-10.166pt}{-2.58514pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}% {rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{30.77252pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{{}{}}}% {{}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{-18.05489pt}{28.19475pt}\pgfsys@invoke{ }\hbox{{\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$0.5$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {{ {{ }}{{{}}} {{ }}{{{}}} {}}} {{ { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } { { { {{ }}{{{}}} {{ }}{{{}}} {}}{}} } {{ { {{ }}{{{}}} {{ }}{{{}}} {}}{ {{}{}} }{}}} } } { } {{}{{}}} {{}{{}}} \pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ % } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ } {}{{}}{} {}{{}}\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{% 0.5,0.5,0.5}\pgfsys@color@gray@stroke{0.5}\pgfsys@invoke{ }\pgfsys@roundcap% \pgfsys@invoke{ }{}\pgfsys@moveto{-3.5pt}{61.55244pt}\pgfsys@stroke% \pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{{}{}}}% {{}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{-10.166pt}{58.97467pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}% {rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$1$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@beginscope\pgfsys@invoke{ }{} \pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} {}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} {}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} \pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} {}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} {}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \pgfsys@beginscope\pgfsys@invoke{ }{} {} {{{ {{ }}{{{}}} {}}}}{{}{ {}{}{}{}{}}{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1% .0}{33.19095pt}{-21.12946pt}\pgfsys@invoke{ }\hbox{{\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\small{$\mathit{x}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\pgfsys@beginscope% \pgfsys@invoke{ }{} {} {{{ {{ }}{{{}}} {}}}}{{ {}{}{}{}{}}{}{}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{% { {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{0.0}{1.0}{-1.0}{% 0.0}{-24.83789pt}{23.06683pt}\pgfsys@invoke{ }\hbox{{\definecolor{% pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }% \pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\small{$\varsigma(x)$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{}{}\hss}% \pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}% \lxSVG@closescope\endpgfpicture}} }\end{array}

(Here $\varsigma(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{1+e^{-x}}$ is the sigmoid function, as illustrated.) We can use these functions to build a network as follows (see also Fig. 2):

\mathrm{comp}\langle\mathrm{layer}_{m}(\mathrm{neuron}_{k}),\mathrm{comp}% \langle\mathrm{layer}_{n}(\mathrm{neuron}_{m}),\mathrm{neuron}_{n}\rangle% \rangle:\boldsymbol{(}\mathbf{real}^{k}\boldsymbol{\mathop{*}}P\boldsymbol{)}% \to\mathbf{real}

(1)

Figure 2. The network in (1) with

k

inputs and two hidden layers.

Here $P\cong\mathbf{real}^{p}$ with $p=(m(k{+}1){+}n(m{+}1){+}n{+}1)$ . This program (1) describes a smooth (infinitely differentiable) function. The goal of automatic differentiation is to find its derivative.

If we $\beta$ -reduce all the $\lambda$ ’s, we end up with a very long function expression just built from the sigmoid function and linear algebra. We can then find a program for calculating its derivative by applying the chain rule. However, automatic differentiation can also be expressed without first $\beta$ -reducing, in a compositional way, by explaining how higher order functions like $(\mathrm{layer})$ and $(\mathrm{comp})$ propagate derivatives. This paper is a semantic analysis of this compositional approach.

The general idea of denotational semantics is to interpret types as spaces and programs as functions between the spaces. In this paper, we propose to use diffeological spaces and smooth functions [Sou80, IZ13] to this end. These satisfy the following three desiderata:

•

$\mathbb{R}$ is a space, and the smooth functions $\mathbb{R}\to\mathbb{R}$ are exactly the functions that are infinitely differentiable;
•

The set of smooth functions $X\to Y$ between spaces again forms a space, so we can interpret function types.
•

The disjoint union of a sequence of spaces again forms a space, and this enables us to interpret variant types and inductive types, e.g. lists of reals form the space $\biguplus_{i=0}^{\infty}\mathbb{R}^{i}$ .

We emphasise that the most standard formulation of differential geometry, using manifolds, does not support spaces of functions. Diffeological spaces seem to us the simplest notion of space that satisfies these conditions, but there are other candidates [BH11, Sta11]. A diffeological space is in particular a set $X$ equipped with a chosen set of curves $C_{X}\subseteq X^{\mathbb{R}}$ and a smooth map $f:X\to Y$ must be such that if $\gamma\in C_{X}$ then $\gamma;f\in C_{Y}$ . This is reminiscent of the method of logical relations.

1.0.2. From smoothness to automatic derivatives at higher types.

Our denotational semantics in diffeological spaces guarantees that all definable functions are smooth. But we need more than just to know that a definable function happens to have a mathematical derivative: we need to be able to find that derivative.

In this paper we focus on forward mode automatic differentiation methods for computing higher derivatives, which are macro translations on syntax (called $\overrightarrow{\mathcal{D}}$ in Section 3). We are able to show that they are correct, using our denotational semantics.

Here there is one subtle point that is central to our development. Although differential geometry provides established derivatives for first order functions (such as $\mathrm{neuron}$ above), there is no canonical notion of derivative for higher order functions (such as $\mathrm{layer}$ and $\mathrm{comp}$ ) in the theory of diffeological spaces (e.g. [CW14]). We propose a new way to resolve this, by interpreting types as triples $(X,X^{\prime},S)$ where, intuitively, $X$ is a space of inhabitants of the type, $X^{\prime}$ is a space serving as a chosen bundle of tangents (or jets, in the case of higher order derivatives) over $X$ , and $S\subseteq X^{\mathbb{R}}\times X^{\prime\mathbb{R}}$ is a binary relation between curves, informally relating curves in $X$ with their tangent (resp. jet) curves in $X^{\prime}$ . This new model gives a denotational semantics for higher order automatic differentiation on a language with higher order functions.

In Section 4 we boil this new approach down to a straightforward and elementary logical relations argument for the correctness of higher order automatic differentiation. The approach is explained in detail in Section 6. We explore some subtleties of non-uniqueness of derivatives of higher order functions in Section 7.

1.0.3. Related work and context.

AD has a long history and has many implementations. AD was perhaps first phrased in a functional setting in [PS08], and there are now a number of teams working on AD in the functional setting (e.g. [WWE⁺19, SFVPJ19, Ell18]), some providing efficient implementations. Although that work does not involve formal semantics, it is inspired by intuitions from differential geometry and category theory.

This paper adds to a very recent body of work on verified automatic differentiation. In the first order setting, there are recent accounts based on denotational semantics in manifolds [FST19, LYRY20] and based on synthetic differential geometry [CGM19], work making a categorical abstraction [CCG⁺20] and work connecting operational semantics with denotational semantics [AP20, Plo18], as well as work focussing on how to correctly differentiate programs that operate on tensors [BML⁺20] and programs that make use of quantum computing [ZHCW20]. Recently there has also been significant progress at higher types. The work of Brunel et al. [BMP20] and Mazza and Pagani [MP21] give formal correctness proofs for reverse-mode derivatives on a linear $\lambda$ -calculus with a particular operational semantics. The work of Barthe et al. [BCLG20] provides a general discussion of some new syntactic logical relations arguments including one very similar to our syntactic proof of Theorem 3. Sherman et al. [SMC20] discuss a differential programming technique that works at higher types, based on exact real arithmetic and relate it to a computable semantics. We understand that the authors of [CGM19] are working on higher types. Vákár [Vák21, VS21, LNV21] phrase and prove correct a reverse mode AD technique on a higher order language based on a similar gluing technique. Vákár [Vák20] extends a standard $\lambda$ -calculus with type recursion, and proves correct a forward-mode AD on such a higher-order language, also using a gluing argument.

The differential $\lambda$ -calculus [ER03] is related to AD, and explicit connections are made in [MO20, Man12]. One difference is that the differential $\lambda$ -calculus allows the addition of terms at all types, and hence vector space models are suitable to interpret all types. This choice would appear peculiar with the variant and inductive types that we consider here, as the dimension of a disjoint union of spaces is only defined locally.

This paper builds on our previous work [HSV20a, Vák20] in which we gave denotational correctness proofs for forward mode AD algorithms for computing first derivatives. Here, we explain how these techniques extend to methods that calculate higher derivatives.

The Faà di Bruno construction has also been investigated [CS11] in the context of Cartesian differential categories.

The idea of directly calculating higher order derivatives by using automatic differentiation methods that work with Taylor approximations (also known as jets in differential geometry) is well-known [GUW00] and it has recently gained renewed interest [Bet18, BJD19]. So far, such “Taylor-mode AD” methods have only been applied to first order functional languages, however. This paper shows how to extend these higher order AD methods to languages with support for higher order functions and algebraic data types.

The two main methods for implementing AD are operator overloading and, the method used in this paper, source code transformation [VMBBL18]. Taylor-mode AD has been seen to be significantly faster than iterated AD in the context of operator overloading [BJD19] in Jax [FJL18]. There are other notable implementations of forward Taylor-mode [BS96, BS97, Kar01, PS07, WGP16]. Some of them are implemented in a functional language [Kar01, PS07]. Taylor-mode implementations use the rich algebraic structure of derivatives to avoid a lot of redundant computations occurring via iterated first order methods and share of a lot of redundant computations. Perhaps the simplest example to see this is with the sin function, whose iterated derivatives only involve sin, cos, and negation. Importantly, most AD tools have the right complexity up to a constant factor, but this constant is quite important in practice and Taylor-mode helps achieve better performance. Another stunning result of a version of Taylor-mode was achieved in [LMG18], where a gain of performance of up to two orders of magnitude was achieved for computing certain Hessian-vector products using Ricci calculus. In essence, the algorithm used is mixed-mode that is derived via jets in [Bet18]. This is further improved in [LMG20]. Tayor-mode can also be useful for ODE solvers and hence will be important for neural differential equations [CRBD18].

Finally, we emphasise that we have chosen the neural network (1) as our running example mainly for its simplicity. Indeed one would typically use reverse-mode AD to train neural networks in practice. There are many other examples of AD outside the neural networks literature: AD is useful whenever derivatives need to be calculated on high dimensional spaces. This includes optimization problems more generally, where the derivative is passed to a gradient descent method (e.g. [RM51, KW⁺52, Qia99, KB14, DHS11, LN89]). Optimization problems involving higher order functions naturally show up in the calculus of variations and its applications in physics, where one typically looks for a function minimizing a certain integral [GSS00]. Other applications of AD are in advanced integration methods, since derivatives play a role in Hamiltonian Monte Carlo [Nea11, HG14] and variational inference [KTR⁺17]. Second order methods for gradient-descent have also been extensively studied. As the basic second order Newton method requires inverting a high dimentional hessian matrix, several alternatives and approximations have been studied. Some of them still require Taylor-like modes of differentiation and require a matrix-vector product where the matrix resembles the hessian or inverse hessian [KK04, Mar10, Ama12].

1.0.4. Summary of contributions.

We have provided a semantic analysis of higher order automatic differentiation. Our syntactic starting point are higher order forward-mode AD macros on a typed higher order language that extend their well-known first order equivalent (e.g. [SFVPJ19, WWE⁺19, HSV20a]). We present these in Section 3 for function types, and in Section 5 we extend them to inductive types and variants. The main contributions of this paper are as follows.

•

We give a denotational semantics for the language in diffeological spaces, showing that every definable expression is smooth (Section 4).
•

We show correctness of the higher order AD macros by a logical relations argument (Th. 3).
•

We give a categorical analysis of this correctness argument with two parts: a universal property satisfied by the macro in terms of syntactic categories, and a new notion of glued space that abstracts the logical relation (Section 6).
•

We then use this analysis to state and prove a correctness argument at all first order types (Th. 8).

Relation to previous work

This paper extends and develops the paper [HSV20a] presented at the 23rd International Conference on Foundations of Software Science and Computation Structure (FoSSaCS 2020). This version includes numerous elaborations, notably the extension of the definition, semantics and correctness of automatic differentiation methods for computing higher order derivatives (introduced in Section 2.2-2.4) and a novel discussion about derivatives of higher-order functions (Section 7).

2. Rudiments of differentiation: how to calculate with dual numbers and Taylor approximations

2.1. First order differentiation: the chain rule and dual numbers.

We will now recall the definition of gradient of a differentiable function, the goal of AD and, and what it means for AD to be correct. Recall that the derivative of a function $f:\mathbb{R}\to\mathbb{R}$ , if it exists, is a function $\nabla f:\mathbb{R}\to\mathbb{R}$ such that for all $a$ , $\nabla f(a)$ is the gradient of $f$ at $a$ in the sense that the function $x\mapsto f(a)+\nabla f(a)\cdot(x-a)$ gives the best linear approximation of $f$ at $a$ . (The gradient $\nabla f(a)$ is often written $\frac{\mathop{}\!\mathrm{d}f(x)}{\mathop{}\!\mathrm{d}x}(a)$ .)

The chain rule for differentiation tells us that we can calculate $\nabla(f;g)(a)=\nabla f(a)\cdot\nabla g(f(a))$ . In that sense, the chain rule tells us how linear approximations to a function transform under post-composition with another function.

To find $\nabla f$ in a compositional way, using the chain rule, two generalizations are reasonable:

•

We need both $f$ and $\nabla f$ when calculating $\nabla(f;g)$ of a composition $f;g$ , using the chain rule, so we are really interested in the pair $(f,\nabla f):\mathbb{R}\to\mathbb{R}\times\mathbb{R}$ ;
•

In building $f$ we will need to consider functions of multiple arguments, such as $+:\mathbb{R}^{2}\to\mathbb{R}$ , and these functions should propagate derivatives.

Thus we are more generally interested in transforming a function $g:\mathbb{R}^{n}\to\mathbb{R}$ into a function $h:(\mathbb{R}\times\mathbb{R})^{n}\to\mathbb{R}\times\mathbb{R}$ in such a way that for any $f_{1}\dots f_{n}:\mathbb{R}\to\mathbb{R}$ ,

(f_{1},\nabla f_{1},\dots,f_{n},\nabla f_{n});h=((f_{1},\dots,f_{n});g,\nabla(% (f_{1},\dots,f_{n});g))\text{.}

(2)

Computing automatically the program representing $h$ , given a program representing $g$ , is the goal of automatic differentiation. An intuition for $h$ is often given in terms of dual numbers. The transformed function operates on pairs of numbers, $(x,x^{\prime})$ , and it is common to think of such a pair as $x+x^{\prime}\epsilon$ for an ‘infinitesimal’ $\epsilon$ . But while this is a helpful intuition, the formalization of infinitesimals can be intricate, and the development in this paper is focussed on the elementary formulation in (2).

A function $h$ satisfying (2) encodes all the partial derivatives of $g$ . For example, if $g\colon\mathbb{R}^{2}\to\mathbb{R}$ , then with $f_{1}(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}x$ and $f_{2}(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}x_{2}$ , by applying (2) to $x_{1}$ we obtain $h(x_{1},1,x_{2},0)\!=\!(g(x_{1},x_{2}),\frac{\partial g(x,x_{2})}{\partial x}(% x_{1}))$ and similarly $h(x_{1},0,x_{2},1)\!=\!(g(x_{1},x_{2}),\frac{\partial g(x_{1},x)}{\partial x}(% x_{2}))$ . And conversely, if $g$ is differentiable in each argument, then a unique $h$ satisfying (2) can be found by taking linear combinations of partial derivatives, for example:

\textstyle h(x_{1},x_{1}^{\prime},x_{2},x_{2}^{\prime})=(g(x_{1},x_{2}),x_{1}^% {\prime}\cdot\frac{\partial g(x,x_{2})}{\partial x}(x_{1})+x_{2}^{\prime}\cdot% \frac{\partial g(x_{1},x)}{\partial x}(x_{2}))\text{.}

(Here, recall that the partial derivative $\frac{\partial g(x,x_{2})}{\partial x}(x_{1})$ is a particular notation for the gradient $\nabla(g(-,x_{2}))(x_{1})$ , i.e. with $x_{2}$ fixed. )

In summary, the idea of differentiation with dual numbers is to transform a differentiable function $g:\mathbb{R}^{n}\to\mathbb{R}$ to a function $h:\mathbb{R}^{2n}\to\mathbb{R}^{2}$ which captures $g$ and all its partial derivatives. We packaged this up in (2) as an invariant which is useful for building derivatives of compound functions $\mathbb{R}\to\mathbb{R}$ in a compositional way. The idea of (first order) forward mode automatic differentiation is to perform this transformation at the source code level.

We say that a macro for AD is correct if, given a semantic model $sem{-}$ , the program $P$ representing $g=\llbracket P\rrbracket$ is transformed by the macro to a program $P^{\prime}$ representing $h=\llbracket P^{\prime}\rrbracket$ . This means in particular that $P^{\prime}$ computes correct partial derivatives of (the function represented by) $P$ .

Smooth functions.

In what follows we will often speak of smooth functions $\mathbb{R}^{k}\to\mathbb{R}$ , which are functions that are continuous and differentiable, such that their derivatives are also continuous and differentiable, and so on.

2.2. Higher order differentiation: the Faà di Bruno formula and Taylor approximations.

We now generalize the above in two directions:

•

We look for the best local approximations to $f$ with polynomials of some order $R$ , generalizing the above use of linear functions ( $R=1$ ).
•

We can work directly with multivariate functions $\mathbb{R}^{k}\to\mathbb{R}$ instead of functions of one variable $\mathbb{R}\to\mathbb{R}$ ( $k=1$ ).

To make this precise, we recall that, given a smooth function $f:\mathbb{R}^{k}\to\mathbb{R}$ and a natural number $R\geq 0$ , the $R$ -th order Taylor approximation of $f$ at ${a\in\mathbb{R}^{k}}$ is defined in terms of the partial derivatives of $f$ :

	$\displaystyle\mathbb{R}^{k}\;\;\quad$	$\displaystyle\to\qquad\mathbb{R}$
	$\displaystyle{x}\qquad$	$\displaystyle\mapsto\sum_{\left\{(\alpha_{1},\ldots,\alpha_{k})\in\mathbb{N}^{% k}\mid\alpha_{1}+\ldots+\alpha_{k}\leq R\right\}}{\frac{1}{\alpha_{1}!\cdot% \ldots\cdot\alpha_{k}!}\frac{\partial^{\alpha_{1}+\ldots+\alpha_{k}}f(x)}{% \partial x_{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}(a)}\cdot(x_{1}-a% _{1})^{\alpha_{1}}\cdot\ldots\cdot(x_{k}-a_{k})^{\alpha_{k}}.$

This is an $R$ -th order polynomial. Similarly to the case of first order derivatives, we can recover the partial derivatives of $f$ up to the $R$ -th order from its Taylor approximation by evaluating the series at basis vectors. See Section 2.3 below for an example.

Recall that the ordering of partial derivatives does not matter for smooth functions (Schwarz/Clairaut’s theorem). So there will be $\binom{R+k-1}{k-1}$ $R$ -th order partial derivatives, and altogether there are $\binom{R+k}{k}$ summands in the $R$ -th order Taylor approximation. (This can be seen by a ‘stars-and-bars’ argument.)

Since there are ${\binom{R+k}{k}}$ partial derivatives of $f_{i}$ of order $\leq R$ , we can store them in the Euclidean space $\mathbb{R}^{\binom{R+k}{k}}$ , which can also be regarded as the space of $k$ -variate polynomials of degree $\leq R$ .

We use a convention of coordinates $\left(y_{\alpha_{1}...\alpha_{k}}\in\mathbb{R}\right)_{(\alpha_{1},\ldots,% \alpha_{k})\in\left\{(\alpha_{1},\ldots,\alpha_{k})\in\mathbb{N}^{k}\mid 0\leq% \alpha_{1}+\ldots+\alpha_{k}\leq R\right\}}$ where $y_{\alpha_{1}\ldots\alpha_{k}}$ is intended to represent a partial derivative $\frac{\partial^{\alpha_{1}+...+\alpha_{k}}f}{\partial x_{1}^{\alpha_{1}}\cdots% \partial x_{k}^{\alpha_{k}}}(a)$ for some function $f:\mathbb{R}^{k}\to\mathbb{R}$ . We will choose these coordinates in lexicographic order of the multi-indices $(\alpha_{1},\ldots,\alpha_{k})$ , that is, the indexes in the Euclidean space $\mathbb{R}^{\binom{R+k}{k}}$ will typically range from $(0,\ldots,0)$ to $(R,0,\ldots,0)$ .

The $(k,R)$ -Taylor representation of a function $g:\mathbb{R}^{n}\to\mathbb{R}$ is a function $h:\left(\mathbb{R}^{\binom{R+k}{k}}\right)^{n}\to\mathbb{R}^{\binom{R+k}{k}}$ that transforms the partial derivatives of $f:\mathbb{R}^{k}\to\mathbb{R}^{n}$ of order $\leq R$ under postcomposition with $g$ :

{\left({\left(\frac{\partial^{\alpha_{1}+\ldots+\alpha_{k}}f_{j}(x)}{\partial x% _{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}\right)}_{(\alpha_{1},...,% \alpha_{k})=(0,...,0)}^{(R,0,...,0)}\right)}_{j=1}^{n};h={\left({\left(\frac{% \partial^{\alpha_{1}+\ldots+\alpha_{k}}((f_{1},\ldots,f_{n});g)(x)}{\partial x% _{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}\right)}_{(\alpha_{1},...,% \alpha_{k})=(0,...,0)}^{(R,0,...,0)}\right)}_{j=1}^{n}\text{.}

(3)

Thus the Taylor representation generalizes the dual numbers representation ( $R=k=1$ ).

To explicitly calculate the Taylor representation for a smooth function, we recall a generalization of the chain rule to higher derivatives. The chain rule tells us how the coefficients of linear approximations transform under composition of the functions. The Faà di Bruno formula [Sav06, EM03, CS96] tells us how coefficients of Taylor approximations – that is, higher derivatives – transform under composition. We recall the multivariate form from [Sav06, Theorem 2.1]. Given functions $f=(f_{1},\ldots,f_{l}):\mathbb{R}^{k}\to\mathbb{R}^{l}$ and $g:\mathbb{R}^{l}\to\mathbb{R}$ , for $\alpha_{1}+\ldots+\alpha_{k}>0$ ,

	$\displaystyle\frac{\partial^{\alpha_{1}+\ldots+\alpha_{k}}(f;g)(x)}{\partial x% _{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}(a)$	$\displaystyle=\alpha_{1}!\cdot\ldots\cdot\alpha_{k}!\cdot\sum_{\left\{(\beta_{% 1},\ldots,\beta_{l})\in\mathbb{N}^{l}\mid 1\leq\beta_{1}+\ldots+\beta_{l}\leq% \alpha_{1}+\ldots+\alpha_{k}\right\}}\frac{\partial^{\beta_{1}+\ldots+\beta_{l% }}g(y)}{\partial y_{1}^{\beta_{1}}\cdots\partial y_{l}^{\beta_{l}}}(f(a))\cdot$
		$\displaystyle\sum_{\left\{((e^{1}_{1},\ldots,e^{1}_{l}),\ldots,(e^{q}_{1},% \ldots,e^{q}_{l}))\in(\mathbb{N}^{l})^{q}\mid e_{j}^{1}+\ldots+e_{j}^{q}=\beta% _{j},(e^{1}_{1}+\ldots+e^{1}_{l})\cdot\alpha^{1}_{i}+\ldots+(e^{q}_{1}+\ldots+% e^{q}_{l})\cdot\alpha^{q}_{i}=\alpha_{i}\right\}}$
		$\displaystyle\prod_{r=1}^{q}\prod_{j=1}^{l}\frac{1}{e^{r}_{j}!}\left(\frac{1}{% \alpha^{r}_{1}!\cdot\ldots\cdot\alpha^{r}_{k}!}\frac{\partial^{\alpha^{r}_{1}+% \cdots+\alpha^{r}_{k}}f_{j}(x)}{\partial^{\alpha^{r}_{1}}x_{1}\cdots\partial^{% \alpha^{r}_{k}}x_{k}}(a)\right)^{e^{r}_{j}},$

where $(\alpha^{1}_{1},\ldots,\alpha^{1}_{k}),\ldots,(\alpha^{q}_{1},\ldots,\alpha^{q% }_{k})\in\mathbb{N}^{k}$ are an enumeration of all the vectors $(\alpha^{r}_{1},\ldots,\alpha^{r}_{k})$ of $k$ natural numbers such that $\alpha^{r}_{j}\leq\alpha_{j}$ and $\alpha^{r}_{1}+\ldots+\alpha^{r}_{k}>0$ and we write $q$ for the number of such vectors. The details of this formula reflect the complicated combinatorics the arise from repeated applications of the chain and product rules for differentiation that one uses to prove it. Conceptually, however, it is rather straightforward: it tells us that the coefficients of the $R$ -th order Taylor approximation of $f;g$ can be expressed exclusively in terms of those of $f$ and $g$ .

Thus the Faà di Bruno formula uniquely determines the Taylor approximation $h:\left(\mathbb{R}^{\binom{R+k}{k}}\right)^{n}\to\mathbb{R}^{\binom{R+k}{k}}$ in terms of the derivatives of $g:\mathbb{R}^{n}\to\mathbb{R}$ of order $\leq R$ , and we can also recover all such derivatives from $h$ .

2.3. Example: a two-dimensional second order Taylor series

As an example, we can specialize the Faà di Bruno formula above to the second order Taylor series of a function $f:\mathbb{R}^{2}\to\mathbb{R}^{l}$ and its behaviour under postcomposition with a smooth function $g:\mathbb{R}^{l}\to\mathbb{R}$ :

\displaystyle\frac{\partial^{2}(f;g)(x)}{\partial x_{i}\partial x_{i^{\prime}}% }(a)

\displaystyle=\sum_{j=1}^{l}\frac{\partial g(y)}{\partial y_{j}}(f(a))\frac{% \partial^{2}f_{j}(x)}{\partial x_{i}\partial x_{i^{\prime}}}(a)+\sum_{j,j^{% \prime}=1}^{l}\frac{\partial^{2}g(y)}{\partial y_{j}\partial y_{j^{\prime}}}(f% (a))\frac{\partial f_{j^{\prime}}(x)}{\partial x_{i}}(a)\frac{\partial f_{j}(x% )}{\partial x_{i^{\prime}}}(a),

where $i,i^{\prime}\in\left\{1,2\right\}$ might either coincide or be distinct.

Rather than working with the full $(2,2)$ -Taylor representation of $g$ , we ignore the non-mixed second order derivatives $y_{02}^{j}=\frac{\partial^{2}f_{j}(x)}{\partial x_{2}^{2}}$ and $y_{20}^{j}=\frac{\partial^{2}f_{j}(x)}{\partial x_{1}^{2}}$ for the moment, and we represent the derivatives of order $\leq 2$ of $f_{j}:\mathbb{R}^{2}\to\mathbb{R}$ (at some point $a$ ) as the numbers

(y_{00}^{j},y_{01}^{j},y_{10}^{j},y_{11}^{j})=\left(f_{j}(a),\frac{\partial f_% {j}(x)}{\partial x_{2}}(a),\frac{\partial f_{j}(x)}{\partial x_{1}}(a),\frac{% \partial^{2}f_{j}(x)}{\partial x_{1}\partial x_{2}}(a)\right)\in\mathbb{R}^{4}

and we can choose a similar representation for the derivatives of $(f;g)$ . Then, we observe that the Faà di Bruno formula induces the function $h:(\mathbb{R}^{4})^{l}\to\mathbb{R}^{4}$

	$\displaystyle h((y_{00}^{1},y_{01}^{1},y_{10}^{1},y_{11}^{1}),\ldots,(y_{00}^{% l},y_{01}^{l},y_{10}^{l},y_{11}^{l}))=$
	$\displaystyle\left(\begin{array}[]{l}g(y_{00}^{1},\ldots,y_{00}^{l})\\ \sum_{j=1}^{l}\frac{\partial g(y^{1},\ldots,y^{l})}{\partial y_{j}}(y_{00}^{1}% ,\ldots,y_{00}^{l})\cdot y_{01}^{j}\\ \sum_{j=1}^{l}\frac{\partial g(y^{1},\ldots,y^{l})}{\partial y_{j}}(y_{00}^{1}% ,\ldots,y_{00}^{l})\cdot y_{10}^{j}\\ \sum_{j=1}^{l}\frac{\partial g(y^{1},\ldots,y^{l})}{\partial y_{j}}(y_{00}^{1}% ,\ldots,y_{00}^{l})\cdot y_{11}^{j}+\sum_{j,j^{\prime}=1}^{l}\frac{\partial^{2% }g(y^{1},\ldots,y^{l})}{\partial y_{j}\partial y_{j^{\prime}}}(y_{00}^{1},% \ldots,y_{00}^{l})\cdot y_{10}^{j}\cdot y_{01}^{j^{\prime}}\end{array}\right).$

In particular, we can note that

\displaystyle h((y_{00}^{1},y_{01}^{1},y_{10}^{1},0),\ldots,(y_{00}^{l},y_{01}% ^{l},y_{10}^{l},0))=\left(\begin{array}[]{l}g(y_{00}^{1},\ldots,y_{00}^{l})\\ \sum_{j=1}^{l}\frac{\partial g(y^{1},\ldots,y^{l})}{\partial y_{j}}(y_{00}^{1}% ,\ldots,y_{00}^{l})\cdot y_{01}^{j}\\ \sum_{j=1}^{l}\frac{\partial g(y^{1},\ldots,y^{l})}{\partial y_{j}}(y_{00}^{1}% ,\ldots,y_{00}^{l})\cdot y_{10}^{j}\\ \sum_{j,j^{\prime}=1}^{l}\frac{\partial^{2}g(y^{1},\ldots,y^{l})}{\partial y_{% j}\partial y_{j^{\prime}}}(y_{00}^{1},\ldots,y_{00}^{l})\cdot y_{10}^{j}\cdot y% _{01}^{j^{\prime}}\end{array}\right).

We see can use this method to calculate any directional first and second order derivative of $g$ in one pass. For example, if $l=3$ , so $g:\mathbb{R}^{3}\to\mathbb{R}$ , then the last component of $h((x,x^{\prime},x^{\prime\prime},0),(y,y^{\prime},y^{\prime\prime},0),(z,z^{% \prime},z^{\prime\prime},0))$ is the result of taking the first derivative in direction $(x^{\prime},y^{\prime},z^{\prime})$ and the second derivative in direction $(x^{\prime\prime},y^{\prime\prime},z^{\prime\prime})$ , and evaluating at $(x,y,z)$ .

In the proper Taylor representation we explicitly include the non-mixed second order derivatives as inputs and outputs, leading to a function $h^{\prime}:(\mathbb{R}^{6})^{l}\to\mathbb{R}^{6}$ . Above we have followed a common trick to avoid some unnecessary storage and computation, since these extra inputs and outputs are not required for computing the second order derivatives of $g$ . For instance, if $l=2$ then the last component of $h((x,1,1,0),(y,0,0,0))$ computes $\frac{\partial^{2}g(x,y)}{\partial x^{2}}(x,y)$ .

2.4. Example: a one-dimensional second order Taylor series

As opposed to (2,2)-AD, (1,2)-AD computes the first and second order derivatives in the same direction. For example, if $g:\mathbb{R}^{2}\to\mathbb{R}$ is a smooth function, then $h:(\mathbb{R}^{3})^{2}\to\mathbb{R}^{3}$ . An intuition for $h$ can be given in terms of triple numbers. The transformed function operates on triples of numbers, $(x,x^{\prime},x^{\prime\prime})$ , and it is common to think of such a triple as $x+x^{\prime}\epsilon+x^{\prime\prime}\epsilon^{2}$ for an ‘infinitesimal’ $\epsilon$ which has the property that $\epsilon^{3}=0$ . For instance we have

	$\displaystyle h((x_{1},1,0),(x_{2},0,0))=(g(x_{1},x_{2}),\frac{\partial g(x,x_% {2})}{\partial x}(x_{1}),\frac{\partial^{2}g(x,x_{2})}{\partial x^{2}}(x_{1}))$
	$\displaystyle h((x_{1},0,0),(x_{2},1,0))=(g(x_{1},x_{2}),\frac{\partial g(x_{1% },x)}{\partial x}(x_{2}),\frac{\partial^{2}g(x_{1},x)}{\partial x^{2}}(x_{2}))$
	$\displaystyle h((x_{1},1,0),(x_{2},1,0))=(g(x_{1},x_{2}),\frac{\partial g(x,x_% {2})}{\partial x}(x_{1})+\frac{\partial g(x_{1},x)}{\partial x}(x_{2}),$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\frac{\partial^{2}% g(x,x_{2})}{\partial x^{2}}(x_{1})+\frac{\partial^{2}g(x_{1},x)}{\partial x^{2% }}(x_{2})+2\frac{\partial^{2}g(x,y)}{\partial x\partial y}(x_{1},x_{2}))$

We see that we directly get non-mixed second-order partial derivatives but not the mixed-ones. We can recover $\frac{\partial^{2}g(x,y)}{\partial x\partial y}(x_{1},x_{2})$ as $\frac{1}{2}(h((x_{1},1,0),(x_{2},1,0))-h((x_{1},1,0),(x_{2},0,0))-h((x_{1},0,0% ),(x_{2},1,0)))$ .

More generally, if $g:\mathbb{R}^{l}\to\mathbb{R}$ , then $h:(\mathbb{R}^{3})^{l}\to\mathbb{R}^{3}$ satisfies:

\displaystyle h((x_{1},x^{\prime}_{1},0),\ldots,(x_{l},x^{\prime}_{l},0))=% \left(\begin{array}[]{l}g(x_{1},\ldots,x_{l})\\ \sum_{i=1}^{l}\frac{\partial g(x_{1},\ldots,x_{l})}{\partial x_{i}}(x_{1},% \ldots,x_{l})\cdot x^{\prime}_{i}\\ \sum_{i,j=1}^{l}\frac{\partial^{2}g(x_{1},\ldots,x_{l})}{\partial x_{i}% \partial x_{j}}(x_{1},\ldots,x_{l})\cdot x^{\prime}_{i}\cdot x^{\prime}_{j}% \end{array}\right).

We can always recover the mixed second order partial derivatives from this but this requires several computations involving $h$ . This is thus different from the (2,2) method which was more direct.

2.5. Remark

In the rest of this article, we study forward-mode $(k,R)$ -automatic differentiation for a language with higher-order functions. The reader may like to fix $k=R=1$ for a standard automatic differentiation with first-order derivatives, based on dual numbers. This is the approach taken in the conference version of this paper [HSV20b]. But the generalization to higher-order derivatives with arbitrary $k$ and $R$ flows straightforwardly through the whole narrative.

3. A Higher Order Forward-Mode AD Translation

3.1. A simple language of smooth functions.

We consider a standard higher order typed language with a first order type $\mathbf{real}$ of real numbers. The types $({\tau},{\sigma})$ and terms $({t},{s})$ are as follows.

$\begin{array}[t]{l@{\quad\!\!}*3{l@{}}@{\,}l}{\tau},{\sigma},{\rho}&::=&&% \mspace{-25.0mu}\qquad\text{types}\\ &\mathrel{\lvert}&\mathbf{real}&\qquad\text{real numbers}\\ &\mathrel{\lvert}&\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}\dots% \boldsymbol{\mathop{*}}{\tau}_{n}\boldsymbol{)}&\qquad\text{finite product}\\ &\mathrel{\lvert}&{\tau}\to{\sigma}&\qquad\text{function}\\[6.0pt] {t},{s},{r}&::=&&\mspace{-25.0mu}\qquad\text{terms}\\ &&{x}&\qquad\text{variable}\\ &\mathrel{\lvert}&\mathsf{op}({t}_{1},\ldots,{t}_{n})&\qquad\text{operations (% including constants)}\\ &\mathrel{\lvert}&\langle{t}_{1},\dots,{t}_{n}\rangle\ \mathrel{\lvert}\mathbf% {case}\,{t}\,\mathbf{of}\,\langle{x}_{1},\dots,{x}_{n}\rangle\to{s}&\qquad% \text{tuples/pattern matching}\\ &\mathrel{\lvert}&\lambda{x}.{t}\ \mathrel{\lvert}{t}\,{s}&\qquad\text{% function abstraction/application}\\ \end{array}$

The typing rules are in Figure 3. We have included some abstract basic $n$ -ary operations $\mathsf{op}\in\mathsf{Op}_{n}$ for every $n\in\mathbb{N}$ . These are intended to include the usual (smooth) mathematical operations that are used in programs to which automatic differentiation is applied. For example,

•

for any real constant $c\in\mathbb{R}$ , we typically include a constant $\underline{c}\in\mathsf{Op}_{0}$ ; we slightly abuse notation and will simply write $\underline{c}$ for $\underline{c}()$ in our examples;
•

we include some unary operations such as $\varsigma\in\mathsf{Op}_{1}$ which we intend to stand for the usual sigmoid function, $\varsigma(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{1+e^{-x}}$ ;
•

we include some binary operations such as addition and multiplication $(+),(*)\in\mathsf{Op}_{2}$ ;

We add some simple syntactic sugar $t-u\stackrel{{\scriptstyle\mathrm{def}}}{{=}}t+\underline{(-1)}*u$ and, for some natural number $n$ ,

n\cdot{t}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\overbrace{{t}+...+{t}}^{% \text{$n$ times}}\qquad\text{and}\qquad{t}^{n}\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\overbrace{{t}*...*{t}}^{\text{$n$ times}}

Similarly, we will frequently denote repeated sums and products using $\sum$ - and $\prod$ -signs, respectively: for example, we write ${t}_{1}+...+{t}_{n}$ as $\sum_{i\in\{1,...,n\}}{t}_{i}$ and ${t}_{1}*...*{t}_{n}$ as $\prod_{i\in\left\{1,...,n\right\}}{t}_{i}$ . This in addition to programming sugar such as $\mathbf{let}\,{x}=\,{t}\,\mathbf{in}\,{s}$ for $(\lambda{x}.{{s}})\,{t}$ and $\lambda\langle{x}_{1},\ldots,{x}_{n}\rangle.{{t}}$ for $\lambda{x}.{\mathbf{case}\,{x}\,\mathbf{of}\,\langle{x}_{1},\ldots,{x}_{n}% \rangle\to{t}}$ .

$\begin{array}[]{c}\inferrule{\Gamma\vdash{t}_{1}:\mathbf{real}\;\;\dots\;\;% \Gamma\vdash{t}_{n}:\mathbf{real}}{\Gamma\vdash\mathsf{op}({t}_{1},\ldots,{t}_% {n}):\mathbf{real}}(\mathsf{op}\in\mathsf{Op}_{n})\\[12.0pt] \\ \inferrule{\Gamma\vdash{t}_{1}:{\tau}_{1}\;\;\dots\;\;\Gamma\vdash{t}_{n}:{% \tau}_{n}}{\Gamma\vdash\langle{t}_{1},\dots,{t}_{n}\rangle:\boldsymbol{(}{\tau% }_{1}\boldsymbol{\mathop{*}}\dots\boldsymbol{\mathop{*}}{\tau}_{n}\boldsymbol{% )}}\qquad\inferrule{\Gamma\vdash{t}:\boldsymbol{(}{\sigma}_{1}\boldsymbol{% \mathop{*}}\dots\boldsymbol{\mathop{*}}{\sigma}_{n}\boldsymbol{)}\;\;\Gamma,{{% x}_{1}\colon{\sigma}_{1},{.}{.}{.},{x}_{n}\colon{\sigma}_{n}}\vdash{s}:{\tau}}% {\Gamma\vdash\mathbf{case}\,{t}\,\mathbf{of}\,\langle{x}_{1},\dots,{x}_{n}% \rangle\to{s}:{\tau}}\\[12.0pt] \\ \inferrule{~{}}{\Gamma\vdash{x}:{\tau}}(({x}:{\tau})\in\Gamma)\qquad\inferrule% {\Gamma,{x}:{\tau}\vdash{t}:{\sigma}}{\Gamma\vdash\lambda{x}:{\tau}.{t}:{\tau}% \to{\sigma}}\qquad\inferrule{\Gamma\vdash{t}:{\sigma}\to{\tau}\\ \Gamma\vdash{s}:{\sigma}}{\Gamma\vdash{t}\,{s}:{\tau}}\end{array}$

Figure 3. Typing rules for the simple language.

3.2. Syntactic automatic differentiation: a functorial macro.

The aim of higher order forward mode AD is to find the $(k,R)$ -Taylor representation of a function by syntactic manipulations, for some choice of $(k,R)$ that we fix. For our simple language, we implement this as the following inductively defined macro $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ on both types and terms (see also [WWE⁺19, SFVPJ19]). For the sake of legibility, we simply write $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ as $\overrightarrow{\mathcal{D}}$ here and leave the dimension $k$ and order $R$ of the Taylor representation implicit. The following definition is for general $k$ and $R$ , but we treat specific cases afterwards in Example 3.2.

	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}\to{\sigma})% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}\scalebox{0.8}{$\overrightarrow{% \mathcal{D}}$}({\tau})\to\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({% \sigma})\qquad\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}_{1}% \boldsymbol{\mathop{}}...\boldsymbol{\mathop{}}{\tau}_{n})\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}% ({\tau}_{1})}\boldsymbol{\mathop{}}...\boldsymbol{\mathop{}}{\scalebox{0.8}{% $\overrightarrow{\mathcal{D}}$}({\tau}_{n})}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\mathbf{real})% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}{\mathbf{real}}^{\binom{R+k}{k}}% \quad\text{(i.e.,~{}the type of tuples of reals of length $\textstyle\binom{R+% k}{k}$)}$

	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({x})\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}{x}\hskip 80.0pt\scalebox{0.8}{$\overrightarrow% {\mathcal{D}}$}(\underline{c})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \langle\underline{c},0\rangle$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\lambda{x}.{t})% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda{x}.{\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({t})}\hskip 12.0pt\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({t}\,{s})\stackrel{{\scriptstyle\mathrm{def}}}{% {=}}\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t})\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({s})\hskip 12.0pt\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}(\langle{t}_{1},\dots,{t}_{n}\rangle)\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}\langle\scalebox{0.8}{$\overrightarrow{\mathcal% {D}}$}({t}_{1}),\dots,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t}_{n})\rangle$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\mathbf{case}\,{t% }\,\mathbf{of}\,\langle{x}_{1},\dots,{x}_{n}\rangle\to{s}})\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{% \mathcal{D}}$}({t})\,\mathbf{of}\,\langle{x}_{1},\dots,{x}_{n}\rangle\to% \scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({s})$
	$\displaystyle\begin{array}[]{ll}\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}% (\mathsf{op}({t}_{1},\ldots,{t}_{n}))\stackrel{{\scriptstyle\mathrm{def}}}{{=}% }{}&\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t}_{1})\,% \mathbf{of}\,\langle{x}_{0...0}^{1},...,{x}_{R,0...0}^{1}\rangle\to\\ &\vdots\\ &{\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t}_{n})\,% \mathbf{of}\,\langle{x}_{0...0}^{n},...,{x}_{R,0...0}^{n}\rangle\to}\\ &\begin{array}[]{lll}\langle&D^{0...0}\mathsf{op}(x_{0...0}^{1},...,x_{R,0...0% }^{1},...,x_{0...0}^{n},....,x_{...0}^{n}),\\ &\cdots,\\ &D^{R...0}\mathsf{op}(x_{0...0}^{1},...,x_{R,0...0}^{1},...,x_{0...0}^{n},....% ,x_{R,0...0}^{n})&\rangle\end{array}\end{array}$

	where	$\displaystyle D^{0...0}\mathsf{op}(x_{0...0}^{1},...,x_{R,0...0}^{1},...,x_{0.% ..0}^{n},....,x_{R,0...0}^{n})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \mathsf{op}(x_{0...0}^{1},...,x_{0...0}^{n})$
		$\displaystyle\begin{array}[]{@{}l}D^{\alpha_{1}....\alpha_{k}}\mathsf{op}(x_{0% ...0}^{1},...,x_{R,0...0}^{1},...,x_{0...0}^{n},....,x_{R,0...0}^{n})\stackrel% {{\scriptstyle\mathrm{def}}}{{=}}~{}\qquad\text{\scriptsize(for $\alpha_{1}+..% .+\alpha_{k}>0$)}\\ \begin{array}[]{c}{\alpha_{1}!\cdot\ldots\cdot\alpha_{k}!}\cdot\sum_{\left\{(% \beta_{1},\ldots,\beta_{n})\in\mathbb{N}^{l}\mid 1\leq\beta_{1}+\ldots+\beta_{% n}\leq\alpha_{1}+\ldots+\alpha_{k}\right\}}\partial_{\beta_{1}\cdots\beta_{n}}% \mathsf{op}({x}_{0...0}^{1},\ldots,{x}_{0...0}^{n})*\\ \sum_{\left\{((e^{1}_{1},\ldots,e^{1}_{l}),\ldots,(e^{q}_{1},\ldots,e^{q}_{l})% )\in(\mathbb{N}^{l})^{q}\mid e_{j}^{1}+\ldots+e_{j}^{q}=\beta_{j},(e^{1}_{1}+% \ldots+e^{1}_{l})\cdot\alpha^{1}_{i}+\ldots+(e^{q}_{1}+\ldots+e^{q}_{l})\cdot% \alpha^{q}_{i}=\alpha_{i}\right\}}\\ \prod_{r=1}^{q}\prod_{j=1}^{l}{\frac{1}{e^{r}_{j}!}}\cdot\left({\frac{1}{% \alpha^{r}_{1}!\cdot\ldots\cdot\alpha^{r}_{k}!}}\cdot{x}_{\alpha_{1}\cdots% \alpha_{k}}^{j}\right)^{e^{r}_{j}}\text{.}\end{array}\end{array}$

Here, $(\partial_{\beta_{1}\cdots\beta_{n}}\mathsf{op})(x_{1},\ldots,x_{n})$ are some chosen terms of type $\mathbf{real}$ in the language with free variables from $x_{1},\ldots,x_{n}$ . We think of these terms as implementing the partial derivative $\frac{\partial^{\beta_{1}+...+\beta_{n}}\llbracket\mathsf{op}\rrbracket(x_{1},% ...,x_{n})}{\partial x_{1}^{\beta_{1}}\cdots\partial x_{n}^{\beta_{n}}}$ of the smooth function $\llbracket\mathsf{op}\rrbracket:\mathbb{R}^{n}\to\mathbb{R}$ that $\mathsf{op}$ implements. For example, we could choose the following representations of derivatives of order $\leq 2$ of our example operations

\begin{array}[]{ll}\partial_{01}(+)(x_{1},x_{2})=\underline{1}&\partial_{02}(+% )(x_{1},x_{2})=\underline{0}\\ \partial_{10}(+)(x_{1},x_{2})=\underline{1}&\partial_{11}(+)(x_{1},x_{2})=% \underline{0}\\ \partial_{20}(+)(x_{1},x_{2})=\underline{0}\\[6.0pt] \partial_{01}(*)(x_{1},x_{2})=x_{1}&\partial_{02}(*)(x_{1},x_{2})=\underline{0% }\\ \partial_{10}(*)(x_{1},x_{2})=x_{2}&\partial_{11}(*)(x_{1},x_{2})=\underline{1% }\\ \partial_{20}(*)(x_{1},x_{2})=\underline{0}\\[6.0pt] \partial_{1}(\varsigma)({x})=\mathbf{let}\,{y}=\,\varsigma({x})\,\mathbf{in}\,% {y}*(\underline{1}-{y}){}&\partial_{2}(\varsigma)({x})=\mathbf{let}\,{y}=\,% \varsigma({x})\,\mathbf{in}\\ &\phantom{\partial_{2}(\varsigma)({x})=\;}\mathbf{let}\,{z}=\,{y}*(\underline{% 1}-{y})\,\mathbf{in}\,{z}*(\underline{1}-\underline{2}*{y})\end{array}

Note that our rules, in particular, imply that $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\underline{c})=\langle% \underline{c},\underline{0},\ldots,\underline{0}\rangle$ .

{exa}

[ $(1,1)$ - and $(2,2)$ -AD] Our choices of partial derivatives of the example operations are sufficient to implement $(k,R)$ -Taylor forward AD with $R\leq 2$ . To be explicit, the distinctive formulas for $(1,1)$ - and $(2,2)$ -AD methods (specializing our abstract definition of $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ above) are

	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}(\mathbf{% real})=\boldsymbol{(}\mathbf{real}\boldsymbol{\mathop{*}}\mathbf{real}% \boldsymbol{)}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}(\mathsf{op% }({t}_{1},\ldots,{t}_{n}))=$
	$\displaystyle\qquad\begin{array}[]{l}\mathbf{case}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}_{(1,1)}({t}_{1})\,\mathbf{of}\,\langle{x}^{1}_{% 0},{x}^{1}_{1}\rangle\to\ldots\to\mathbf{case}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}_{(1,1)}({t}_{n})\,\mathbf{of}\,\langle{x}_{0}^{% n},{x}_{1}^{n}\rangle\to\\ \langle\mathsf{op}({x}_{0}^{1},\ldots,{x}_{0}^{n}),\sum_{i=1}^{n}{x}^{i}_{1}*% \partial_{i}\mathsf{op}({x}_{0}^{1},\ldots,{x}_{0}^{n})\rangle\end{array}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)}(\mathbf{% real})=\mathbf{real}^{6}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)}(\mathsf{op% }({t}_{1},...,{t}_{n}))=$
	$\displaystyle\qquad\begin{array}[]{l}\mathbf{case}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}_{(2,2)}({t}_{1})\,\mathbf{of}\,\langle{x}^{1}_{% 00},{x}^{1}_{01},{x}^{1}_{02},{x}^{1}_{10},{x}^{1}_{11},{x}^{1}_{20}\rangle\to% \\ \vdots\\ \mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)}({t}_{n})% \,\mathbf{of}\,\langle{x}^{n}_{00},{x}^{n}_{01},{x}^{n}_{02},{x}^{n}_{10},{x}^% {n}_{11},{x}^{n}_{20}\rangle\to\\ \qquad\begin{array}[]{l}\langle\mathsf{op}({x}^{1}_{00},\ldots,{x}^{n}_{00}),% \\ \sum_{i=1}^{n}{x}^{i}_{01}\partial_{\hat{i}}\mathsf{op}({x}_{00}^{1},\ldots,{% x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{02}\partial_{{\hat{i}}}\mathsf{op}({x}_{00}^{1},\ldots% ,{x}_{00}^{n})+\sum_{i,j=1}^{n}{x}^{i}_{01}{x}^{j}_{01}\partial_{\widehat{{i% ,j}}}\mathsf{op}({x}_{00}^{1},\ldots,{x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{10}\partial_{\hat{i}}\mathsf{op}({x}_{00}^{1},\ldots,{% x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{11}\partial_{{\hat{i}}}\mathsf{op}({x}_{00}^{1},\ldots% ,{x}_{00}^{n})+\sum_{i,j=1}^{n}{x}^{i}_{10}{x}^{j}_{01}\partial_{\widehat{{i% ,j}}}\mathsf{op}({x}_{00}^{1},\ldots,{x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{20}\partial_{{\hat{i}}}\mathsf{op}({x}_{00}^{1},\ldots% ,{x}_{00}^{n})+\sum_{i,j=1}^{n}{x}^{i}_{10}{x}^{j}_{10}*\partial_{\widehat{{i% ,j}}}\mathsf{op}({x}_{00}^{1},\ldots,{x}_{00}^{n})\rangle\end{array}\end{array}$

where we informally write $\hat{i}$ for the one-hot encoding of $i$ (the sequence of length $n$ consisting exclusively of zeros except in position $i$ where it has a $1$ ) and $\widehat{i,j}$ for the two-hot encoding of $i$ and $j$ (the sequence of length $n$ consisting exclusively of zeros except in positions $i$ and $j$ where it has a $1$ if $i\neq j$ and a $2$ if $i=j$ )

As noted in Section 2, it is often unnecessary to include all components of the $(2,2)$ -algorithm, for example when computing a second order directional derivative. In that case, we may define a restricted $(2,2)$ -AD algorithm that drops the non-mixed second order derivatives from the definitions above and defines $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)^{\prime}}(\mathbf{real})% =\mathbf{real}^{4}$ and
$\begin{array}[]{ll}\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)^{% \prime}}(\mathsf{op}({t}_{1},...,{t}_{n}))=&\mathbf{case}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}_{(2,2)^{\prime}}({t}_{1})\,\mathbf{of}\,\langle% {x}^{1}_{00},{x}^{1}_{01},{x}^{1}_{10},{x}^{1}_{11}\rangle\to\\ &\vdots\\ &\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)^{\prime}% }({t}_{n})\,\mathbf{of}\,\langle{x}^{n}_{00},{x}^{n}_{01},{x}^{n}_{10},{x}^{n}% _{11}\rangle\to\\ &\begin{array}[]{l}\langle\mathsf{op}({x}^{1}_{00},\ldots,{x}^{n}_{00}),\\ \sum_{i=1}^{n}{x}^{i}_{01}*\partial_{\hat{i}}\mathsf{op}({x}_{00}^{1},\ldots,{% x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{10}*\partial_{\hat{i}}\mathsf{op}({x}_{00}^{1},\ldots,{% x}_{00}^{n}),\\ \sum_{i=1}^{n}{x}^{i}_{11}*\partial_{{\hat{i}}}\mathsf{op}({x}_{00}^{1},\ldots% ,{x}_{00}^{n})\\ {}+\sum_{i,i^{\prime}=1}^{n}{x}^{i}_{10}*{x}^{i}_{01}*\partial_{\hat{\hat{i}}}% \mathsf{op}({x}_{00}^{1},\ldots,{x}_{00}^{n})\rangle.\end{array}\end{array}$

We extend $\overrightarrow{\mathcal{D}}$ to contexts: $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\{{x}_{1}{:}{\tau}_{1},{.}{.}{.% },{x}_{n}{:}{\tau}_{n}\})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\{{x}_{1}{:% }\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}_{1}),{.}{.}{.},{x}_{n}{% :}\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}_{n})\}$ . This turns $\overrightarrow{\mathcal{D}}$ into a well-typed, functorial macro in the following sense.

Lemma 1 (Functorial macro).

If $\Gamma\vdash{t}:{\tau}$ then $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\Gamma)\vdash\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({t}):\scalebox{0.8}{$\overrightarrow{\mathcal{D% }}$}({\tau})$ .
If $\Gamma,{x}:{\sigma}\vdash{t}:{\tau}$ and $\Gamma\vdash{s}:{\sigma}$ then $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\Gamma)\vdash\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({t}{}[^{{s}}\!/\!_{{x}}])=\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({t}){}[^{\scalebox{0.8}{$\overrightarrow{% \mathcal{D}}$}({s})}\!/\!_{{x}}]$ .

Proof 3.1.

By induction on the structure of typing derviations.

{exa}

[Inner products] Let us write ${\tau}^{n}$ for the $n$ -fold product $\boldsymbol{(}{\tau}\boldsymbol{\mathop{*}}\dots\boldsymbol{\mathop{*}}{\tau}% \boldsymbol{)}$ . Then, given $\Gamma\vdash{t},{s}:\mathbf{real}^{n}$ we can define their inner product

\begin{array}[]{ll}\Gamma\vdash{t}\cdot_{n}{s}\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}&\mathbf{case}\,{t}\,\mathbf{of}\,\langle{z}_{1},\ldots,{z}_{n}% \rangle\to\\ &\mathbf{case}\,{s}\,\mathbf{of}\,\langle{y}_{1},\ldots,{y}_{n}\rangle\to{z}_{% 1}*{y}_{1}+\dots+{z}_{n}*{y}_{n}:\mathbf{real}\end{array}

To illustrate the calculation of $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}$ , let us expand (and $\beta$ -reduce) $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}({t}\cdot_{2}{s})$ :

		$\displaystyle\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1% ,1)}({t})\,\mathbf{of}\,\langle{z}_{1},{z}_{2}\rangle\to\mathbf{case}\,% \scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}({s})\,\mathbf{of}\,% \langle{y}_{1},{y}_{2}\rangle\to$
		$\displaystyle\mathbf{case}\,{z}_{1}\,\mathbf{of}\,\langle{z}_{1,1},{z}_{1,2}% \rangle\to\mathbf{case}\,{y}_{1}\,\mathbf{of}\,\langle{y}_{1,1},{y}_{1,2}\rangle\to$
		$\displaystyle\mathbf{case}\,{z}_{2}\,\mathbf{of}\,\langle{z}_{2,1},{z}_{2,2}% \rangle\to\mathbf{case}\,{y}_{2}\,\mathbf{of}\,\langle{y}_{2,1},{y}_{2,2}\rangle\to$
		$\displaystyle\qquad\langle{z}_{1,1}{y}_{1,1}+{z}_{2,1}{y}_{2,1}\ ,\ {z}_{1,1% }{y}_{1,2}+{z}_{1,2}{y}_{1,1}+{z}_{2,1}{y}_{2,2}+{z}_{2,2}{y}_{2,1}\rangle$
Let us also expand the calculation of $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)^{\prime}}({t}\cdot_{2}{s})$ :
		$\displaystyle\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2% ,2)^{\prime}}({t})\,\mathbf{of}\,\langle{z}_{1},{z}_{2}\rangle\to\mathbf{case}% \,\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(2,2)^{\prime}}({s})\,% \mathbf{of}\,\langle{y}_{1},{y}_{2}\rangle\to$
		$\displaystyle\mathbf{case}\,{z}_{1}\,\mathbf{of}\,\langle{z}_{1},{z}_{1,1}^{% \prime},{z}_{1,2}^{\prime},{z}_{1}^{\prime\prime}\rangle\to\mathbf{case}\,{y}_% {1}\,\mathbf{of}\,\langle{y}_{1},{y}_{1,1}^{\prime},{y}_{1,2}^{\prime},{y}_{1}% ^{\prime\prime},\rangle\to$
		$\displaystyle\mathbf{case}\,{z}_{2}\,\mathbf{of}\,\langle{z}_{2},{z}_{2,1}^{% \prime},{z}_{2,2}^{\prime},{z}_{2}^{\prime\prime},\rangle\to\mathbf{case}\,{y}% _{2}\,\mathbf{of}\,\langle{y}_{2},{y}_{2,1}^{\prime},{y}_{2,2}^{\prime},{y}_{2% }^{\prime\prime},\rangle\to$
		$\displaystyle\langle{z}_{1}{y}_{1}+{z}_{2}{y}_{2},$
		$\displaystyle\qquad\ {z}_{1}{y}_{1,1}^{\prime}+{z}_{1,1}^{\prime}{y}_{1}+{z}% _{2}{y}_{2,1}^{\prime}+{z}_{2,1}^{\prime}{y}_{2},$
		$\displaystyle\qquad\ {z}_{1}{y}_{1,2}^{\prime}+{z}_{1,2}^{\prime}{y}_{1}+{z}% _{2}{y}_{2,2}^{\prime}+{z}_{2,2}^{\prime}{y}_{2},$
		$\displaystyle\qquad\ {z}_{1}^{\prime\prime}{y}_{1}+{z}_{2}^{\prime\prime}{y}% _{2}+{y}_{1}^{\prime\prime}{z}_{1}+{y}_{2}^{\prime\prime}{z}_{2}+$
		$\displaystyle\qquad{y}_{1,1}^{\prime}{y}_{1,2}^{\prime}+{y}_{2,1}^{\prime}{y% }_{2,2}^{\prime}+{z}_{1,1}^{\prime}{z}_{1,2}^{\prime}+{z}_{2,1}^{\prime}{z}_% {2,2}^{\prime}\rangle$

{exa}

[Neural networks] In our introduction, we provided a program (1) in our language to build a neural network out of expressions $\mathrm{neuron},\mathrm{layer},\mathrm{comp}$ ; this program makes use of the inner product of Ex. 3.1. We can similarly calculate the derivatives of deep neural nets by mechanically applying the macro $\overrightarrow{\mathcal{D}}$ .

4. Semantics of differentiation

Consider for a moment the first order fragment of the language in Section 3, with only one type, $\mathbf{real}$ , and no $\lambda$ ’s or pairs. This has a simple semantics in the category of cartesian spaces and smooth maps. Indeed, a term ${x}_{1}\dots x_{n}:\mathbf{real}\vdash{t}:\mathbf{real}$ has a natural reading as a function $\llbracket{t}\rrbracket:\mathbb{R}^{n}\to\mathbb{R}$ by interpreting our operation symbols by the well-known operations on $\mathbb{R}^{n}\to\mathbb{R}$ with the corresponding name. In fact, the functions that are definable in this first order fragment are smooth. Let us write $\mathbf{CartSp}$ for this category of cartesian spaces ( $\mathbb{R}^{n}$ for some $n$ ) and smooth functions.

The category $\mathbf{CartSp}$ has cartesian products, and so we can also interpret product types, tupling and pattern matching, giving us a useful syntax for constructing functions into and out of products of $\mathbb{R}$ . For example, the interpretation of $(\mathrm{neuron}_{n})$ in (1) becomes

\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\xrightarrow{\llbracket\cdot% _{n}\rrbracket\times{\rm id}_{\mathbb{R}}}\mathbb{R}\times\mathbb{R}% \xrightarrow{\llbracket+\rrbracket}\mathbb{R}\xrightarrow{\llbracket\varsigma% \rrbracket}\mathbb{R}

where $\llbracket\cdot_{n}\rrbracket$ , $\llbracket+\rrbracket$ and $\llbracket\varsigma\rrbracket$ are the usual inner product, addition and the sigmoid function on $\mathbb{R}$ , respectively.

Inside this category, we can straightforwardly study the first order language without $\lambda$ ’s, and automatic differentiation. In fact, we can prove the following by plain induction on the syntax:
The interpretation of the (syntactic) forward AD $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t})$ of a first order term ${t}$ equals the usual (semantic) derivative of the interpretation of ${t}$ as a smooth function.

However, as is well-known, the category $\mathbf{CartSp}$ does not support function spaces. To see this, notice that we have polynomial terms

{x}_{1},\ldots,{x}_{d}:\mathbf{real}\vdash\lambda{y}.\,\textstyle\sum_{n=1}^{d% }{x}_{n}{y}^{n}:\mathbf{real}\to\mathbf{real}

for each $d$ , and so if we could interpret $(\mathbf{real}\to\mathbf{real})$ as a Euclidean space $\mathbb{R}^{p}$ then, by interpreting these polynomial expressions, we would be able to find continuous injections $\mathbb{R}^{d}\to\mathbb{R}^{p}$ for every $d$ , which is topologically impossible for any $p$ , for example as a consequence of the Borsuk-Ulam theorem (see Appx. A).

This lack of function spaces means that we cannot interpret the functions $(\mathrm{layer})$ and $(\mathrm{comp})$ from (1) in $\mathbf{CartSp}$ , as they are higher order functions, even though they are very useful and innocent building blocks for differential programming! Clearly, we could define neural nets such as (1) directly as smooth functions without any higher order subcomponents, though that would quickly become cumbersome for deep networks. A problematic consequence of the lack of a semantics for higher order differential programs is that we have no obvious way of establishing compositional semantic correctness of $\overrightarrow{\mathcal{D}}$ for the given implementation of (1).

We now show that every definable function is smooth, and then in Section 4.2 we show that the $\overrightarrow{\mathcal{D}}$ macro witnesses its derivatives.

4.1. Smoothness at higher types and diffeologies

The aim of this section is to introduce diffeological spaces as a semantic model for the simple language in Section 3. By way of motivation, we begin with a standard set theoretic semantics, where types are interpreted as follows

\textstyle\llfloor\mathbf{real}\rrceil\stackrel{{\scriptstyle\mathrm{def}}}{{=% }}\mathbb{R}\qquad\llfloor\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{*}}\dots% \boldsymbol{\mathop{*}}{\tau}_{n}\boldsymbol{)}\rrceil\stackrel{{\scriptstyle% \mathrm{def}}}{{=}}\prod_{i=1}^{n}\llfloor{\tau}_{i}\rrceil\qquad\llfloor\tau% \to\sigma\rrceil\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\llfloor\tau\rrceil% \to\llfloor\sigma\rrceil)

and a term ${x_{1}:{\tau}_{1},\dots,x_{n}:{\tau}_{n}}\vdash t:{\sigma}$ is interpreted as a function $\prod_{i=1}^{n}\llfloor{\tau}_{i}\rrceil\to\llfloor{\sigma}\rrceil$ , mapping a valuation of the context to a result.

We can show that the interpretation of a term $x_{1}:\mathbf{real},\dots,x_{n}:\mathbf{real}\vdash t:\mathbf{real}$ is always a smooth function $\mathbb{R}^{n}\to\mathbb{R}$ , even if it has higher order subterms. We begin with a fairly standard logical relations proof of this, and then move from this to the semantic model of diffeological spaces.

Proposition 2.

If $x_{1}:\mathbf{real},\dots,x_{n}:\mathbf{real}\vdash t:\mathbf{real}$ then the function $\llfloor t\rrceil:\mathbb{R}^{n}\to\mathbb{R}$ is smooth.

Proof 4.1.

For each type ${\tau}$ define a set $Q_{{\tau}}\subseteq[\mathbb{R}^{k}\to\llfloor{\tau}\rrceil]$ by induction on the structure of types:

	$\displaystyle Q_{\mathbf{real}}$	$\displaystyle=\{f:\mathbb{R}^{k}\to\mathbb{R}~{}\|~{}f\text{ is smooth}\}$
	$\displaystyle Q_{\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{}}\dots% \boldsymbol{\mathop{}}{\tau}_{n}\boldsymbol{)}}$	$\displaystyle=\{\textstyle f:\mathbb{R}^{k}\to\prod_{i=1}^{n}\llfloor{\tau}_{i% }\rrceil~{}\|~{}\forall\vec{r}\in\mathbb{R}^{k}.\ \forall i.\ f_{i}(\vec{r})\in Q% _{{\tau}_{i}}\}$
	$\displaystyle Q_{{{\tau}}\to{{\sigma}}}$	$\displaystyle=\{f:\mathbb{R}^{k}\to\llfloor{\tau}\rrceil\to\llfloor{\sigma}% \rrceil~{}\|~{}\forall g\in Q_{{\tau}}.\,\lambda(\vec{r}).\,f(\vec{r})(g(\vec{r% }))\in Q_{{\sigma}}\}$

Now we show the fundamental lemma: if ${x_{1}:{\tau}_{1},\dots,x_{n}:{\tau}_{n}}\vdash u:{\sigma}$ and $g_{1}\in Q_{{\tau}_{1}}\dots g_{n}\in Q_{{\tau}_{n}}$ then $((g_{1}\dots g_{n});\llfloor u\rrceil)\in Q_{{\sigma}}$ . This is shown by induction on the structure of typing derivations. The only interesting step here is that the basic operations ( $+$ , $*$ , $\varsigma$ etc.) are smooth. We deduce the statement of the theorem by putting $u=t$ , $k=n$ , and letting $g_{i}:\mathbb{R}^{n}\to\mathbb{R}$ be the projections.

At higher types, the logical relations $Q$ show that we can only define functions that send smooth functions to smooth functions, meaning that we can never use them to build first order functions that are not smooth. For example, $(\mathrm{comp})$ in (1) has this property.

This logical relations proof suggests to build a semantic model by interpreting types as sets with structure: for each type we have a set $X$ together with a set $Q^{\mathbb{R}^{k}}_{X}\subseteq[\mathbb{R}^{k}\to X]$ of plots. {defi} A diffeological space $(X,\mathcal{P}_{X})$ consists of a set $X$ together with, for each $n$ and each open subset $U$ of $\mathbb{R}^{n}$ , a set $\mathcal{P}_{X}^{U}\subseteq[U\to X]$ of functions, called plots, such that

•

all constant functions are plots;
•

if $f:V\to U$ is a smooth function and $p\in\mathcal{P}_{X}^{U}$ , then $f;p\in\mathcal{P}_{X}^{V}$ ;
•

if $\left(p_{i}\in\mathcal{P}_{X}^{U_{i}}\right)_{i\in I}$ is a compatible family of plots $(x\in U_{i}\cap U_{j}\Rightarrow p_{i}(x)=p_{j}(x))$ and $\left(U_{i}\right)_{i\in I}$ covers $U$ , then the gluing $p:U\to X:x\in U_{i}\mapsto p_{i}(x)$ is a plot.

We call a function $f:X\to Y$ between diffeological spaces smooth if, for all plots $p\in\mathcal{P}_{X}^{U}$ , we have that $p;f\in\mathcal{P}_{Y}^{U}$ . We write $\mathbf{Diff}(X,Y)$ for the set of smooth maps from $X$ to $Y$ . Smooth functions compose, and so we have a category $\mathbf{Diff}$ of diffeological spaces and smooth functions.

A diffeological space is thus a set equipped with structure. Many constructions of sets carry over straightforwardly to diffeological spaces.

{exa}

[Cartesian diffeologies] Each open subset $U$ of $\mathbb{R}^{n}$ can be given the structure of a diffeological space by taking all the smooth functions $V\to U$ as $\mathcal{P}_{U}^{V}$ . Smooth functions from $V\to U$ in the traditional sense coincide with smooth functions in the sense of diffeological spaces [IZ13]. Thus diffeological spaces have a profound relationship with ordinary calculus.

In categorical terms, this gives a full embedding of $\mathbf{CartSp}$ in $\mathbf{Diff}$ . {exa}[Product diffeologies] Given a family $\left(X_{i}\right)_{i\in I}$ of diffeological spaces, we can equip the product $\prod_{i\in I}X_{i}$ of sets with the product diffeology in which $U$ -plots are precisely the functions of the form $\left(p_{i}\right)_{i\in I}$ for $p_{i}\in\mathcal{P}_{X_{i}}^{U}$ .

This gives us the categorical product in $\mathbf{Diff}$ . {exa}[Functional diffeology] We can equip the set $\mathbf{Diff}(X,Y)$ of smooth functions between diffeological spaces with the functional diffeology in which $U$ -plots consist of functions $f:U\to\mathbf{Diff}(X,Y)$ such that $(u,x)\mapsto f(u)(x)$ is an element of $\mathbf{Diff}(U\times X,Y)$ .

This specifies the categorical function object in $\mathbf{Diff}$ .

We can now give a denotational semantics to our language from Section 3 in the category of diffeological spaces. We interpret each type ${\tau}$ as a set $\llbracket{\tau}\rrbracket$ equipped with the relevant diffeology, by induction on the structure of types:

	$\displaystyle\llbracket\mathbf{real}\rrbracket\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\mathbb{R}\qquad\text{with the standard diffeology}$
	$\displaystyle\llbracket\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{}}\dots% \boldsymbol{\mathop{}}{\tau}_{n}\boldsymbol{)}\rrbracket\ \stackrel{{% \scriptstyle\mathrm{def}}}{{=}}\ \textstyle\prod_{i=1}^{n}\llbracket{\tau}_{i}% \rrbracket\quad\text{with the product diffeology}$
	$\displaystyle\llbracket{\tau}\to{\sigma}\rrbracket\stackrel{{\scriptstyle% \mathrm{def}}}{{=}}\mathbf{Diff}(\llbracket{\tau}\rrbracket,\llbracket{\sigma}% \rrbracket)\quad\text{with the functional diffeology}$

A context $\Gamma=({x}_{1}\colon{\tau}_{1}\dots x_{n}\colon{\tau}_{n})$ is interpreted as a diffeological space $\llbracket\Gamma\rrbracket\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\prod_{i=1% }^{n}\llbracket{\tau}_{i}\rrbracket$ . Now well typed terms $\Gamma\vdash{t}:{\tau}$ are interpreted as smooth functions $\llbracket{t}\rrbracket:\llbracket\Gamma\rrbracket\to\llbracket{\tau}\rrbracket$ , giving a meaning for ${t}$ for every valuation of the context. This is routinely defined by induction on the structure of typing derivations once we choose a smooth function $\llbracket\mathsf{op}\rrbracket:\mathbb{R}^{n}\to\mathbb{R}$ to interpret each $n$ -ary operation $\mathsf{op}\in\mathsf{Op}_{n}$ . For example, constants $\underline{c}:\mathbf{real}$ are interpreted as constant functions; and the first order operations ( $+,*,\varsigma$ ) are interpreted by composing with the corresponding functions, which are smooth: e.g., $\llbracket\varsigma(t)\rrbracket(\rho)\stackrel{{\scriptstyle\mathrm{def}}}{{=% }}\varsigma(\llbracket t\rrbracket(\rho))$ , where $\rho\in\llbracket\Gamma\rrbracket$ . Variables are interpreted as $\llbracket{x}_{i}\rrbracket(\rho)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \rho_{i}$ . The remaining constructs are interpreted as follows, and it is straightforward to show that smoothness is preserved.

	$\displaystyle\llbracket\langle{t}_{1},\dots,{t}_{n}\rangle\rrbracket(\rho)% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\llbracket{t}_{1}\rrbracket(\rho),% \dots,\llbracket{t}_{n}\rrbracket(\rho))$	$\displaystyle\llbracket\lambda{x}{:}{\tau}.{{t}}\rrbracket(\rho)(a)\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}\llbracket{t}\rrbracket(\rho,a)\ \text{($a\in% \llbracket{\tau}\rrbracket$)}$
	$\displaystyle\llbracket\mathbf{case}\,{t}\,\mathbf{of}\,\langle{.}{.}{.}% \rangle\to{s}\rrbracket(\rho)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \llbracket{s}\rrbracket(\rho,\llbracket{t}\rrbracket(\rho))$	$\displaystyle\llbracket{t}\,{s}\rrbracket(\rho)\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\llbracket{t}\rrbracket(\rho)(\llbracket{s}\rrbracket(\rho))$

The logical relations proof of Proposition 2 is reminiscent of diffeological spaces. We now briefly remark on the suitability of the axioms of diffeological spaces (Def 4.1) for a semantic model of smooth programs. The first axiom says that we only consider reflexive logical relations. From the perspective of the interpretation, it recognizes in particular that the semantics of an expression of type $(\mathbf{real}\to\mathbf{real})\to\mathbf{real}$ is defined by its value on smooth functions rather than arbitrary arguments. That is to say, the set-theoretic semantics at the beginning of this section, $\llfloor(\mathbf{real}\to\mathbf{real})\to\mathbf{real}\rrceil$ , is different to the diffeological semantics, $\llbracket(\mathbf{real}\to\mathbf{real})\to\mathbf{real}\rrbracket$ . The second axiom for diffeological spaces ensures that the smooth maps in $\mathbf{Diff}(U,X)$ are exactly the plots in $\mathcal{P}_{X}^{U}$ . The third axiom ensures that categories of manifolds fully embed into $\mathbf{Diff}$ ; it will not play a visible role in this paper — in fact, [BCLG20] prove similar results for a simple language like ours by using plain logical relations (over $\mathbf{Set}$ ) and without demanding the diffeology axioms. However, we expect the third axiom to be crucial for programming with other smooth structures or partiality.

4.2. Correctness of AD

We have shown that a term ${x}_{1}\colon\mathbf{real},\dots,{x}_{n}\colon\mathbf{real}\vdash{t}:\mathbf{real}$ is interpreted as a smooth function $\llbracket{t}\rrbracket:\mathbb{R}^{n}\to\mathbb{R}$ , even if $t$ involves higher order functions (like (1)). Moreover, the macro differentiation $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})$ is a function $\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})% \rrbracket:(\mathbb{R}^{\binom{R+k}{k}})^{n}\to\mathbb{R}^{\binom{R+k}{k}}$ (Proposition 1). This enables us to state a limited version of our main correctness theorem:

Theorem 3 (Semantic correctness of $\overrightarrow{\mathcal{D}}$ (limited)).

For any term
${x}_{1}\colon\mathbf{real},\dots,{x}_{n}\colon\mathbf{real}\vdash{t}:\mathbf{real}$ , the function $\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})\rrbracket$ is the $(k,R)$ -Taylor representation (3) of $\llbracket{t}\rrbracket$ . In detail: for any smooth functions $f_{1}\dots f_{n}:\mathbb{R}^{k}\to\mathbb{R}$ ,

	$\displaystyle{\left({\left(\frac{\partial^{\alpha_{1}+\ldots+\alpha_{k}}f_{j}(% x)}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}\right)}_{(% \alpha_{1},...,\alpha_{k})=(0,...,0)}^{(R,0,...,0)}\right)}_{j=1}^{n};% \llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})\rrbracket=$
	$\displaystyle\qquad\quad{\left({\left(\frac{\partial^{\alpha_{1}+\ldots+\alpha% _{k}}((f_{1},\ldots,f_{n});\llbracket{t}\rrbracket)(x)}{\partial x_{1}^{\alpha% _{1}}\cdots\partial x_{k}^{\alpha_{k}}}\right)}_{(\alpha_{1},...,\alpha_{k})=(% 0,...,0)}^{(R,0,...,0)}\right)}_{j=1}^{n}\text{.}$

For instance, if $n=2$ , then $\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}({t})% \rrbracket({x}_{1},1,{x}_{2},0)=\big{(}\llbracket{t}\rrbracket({x}_{1},{x}_{2}% ),\frac{\partial\llbracket{t}\rrbracket({x},{x}_{2})}{\partial{x}}({x}_{1})% \big{)}$ .

Proof 4.2.

We prove this by logical relations. A categorical version of this proof is in Section 6.2.

For each type ${\tau}$ , we define a binary relation $S_{{\tau}}$ between (open) $k$ -dimensional plots in $\llbracket{\tau}\rrbracket$ and (open) $k$ -dimensional plots in $\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({\tau})\rrbracket$ , i.e. $S_{{\tau}}\subseteq\mathcal{P}_{\llbracket{\tau}\rrbracket}^{\mathbb{R}^{k}}% \times\mathcal{P}_{\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(% k,R)}({\tau})\rrbracket}^{\mathbb{R}^{k}}$ , by induction on ${\tau}$ :

$\displaystyle S_{\mathbf{real}}$	$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}$	$\displaystyle\left\{\left(f,{\left(\frac{\partial^{\alpha_{1}+\ldots+\alpha_{k% }}f(x)}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}\right)}_% {(\alpha_{1},...,\alpha_{k})=(0,...,0)}^{(R,0,...,0)}\right)~{}\Big{\|}~{}f:% \mathbb{R}^{k}\to\mathbb{R}\text{ smooth}\right\}$
$\displaystyle S_{\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{}}...\boldsymbol% {\mathop{}}{\tau}_{n}\boldsymbol{)}}$	$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}$	$\displaystyle\{((f_{1},...,f_{n}),(g_{1},...,g_{n}))\mid(f_{1},g_{1})\in S_{{% \tau}_{1}},...,(f_{n},g_{n})\in S_{{\tau}_{n}}\}$
$\displaystyle S_{{\tau}\to{\sigma}}$	$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}$	$\displaystyle\{(f_{1},f_{2})\mid\forall(g_{1},g_{2})\in S_{{\tau}}.(x{\mapsto}% f_{1}(x)(g_{1}(x)),x{\mapsto}f_{2}(x)(g_{2}(x)))\in S_{{\sigma}}\}$

Then, we establish the following ‘fundamental lemma’:

If ${x}_{1}{:}{\tau}_{1},{.}{.}{.},{x}_{n}{:}{\tau}_{n}\vdash{t}:{\sigma}$ and, for all $1{\leq}i{\leq}n$ , $f_{i}:\mathbb{R}^{k}\to\llbracket{\tau}_{i}\rrbracket$ and
$g_{i}:\mathbb{R}^{k}\to~{}\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D% }}$}_{(k,R)}({\tau}_{i})\rrbracket$ are such that $(f_{i},g_{i})$ is in $S_{{\tau}_{i}}$ , then we have that

$\Big{(}(f_{1},\ldots,f_{n});\llbracket{t}\rrbracket,(g_{1},\ldots,g_{n});% \llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})% \rrbracket\Big{)}$

is in $S_{{\sigma}}$ .

This is proved routinely by induction on the typing derivation of ${t}$ . The case for $\mathsf{op}({t}_{1},\ldots,{t}_{n})$ relies on the precise definition of $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}(\mathsf{op}({t}_{1},% \ldots,{t}_{n}))$ .

We conclude the theorem from the fundamental lemma by considering the case where ${\tau}_{i}={\sigma}=\mathbf{real}$ , $m=n$ and $s_{i}=y_{i}$ .

5. Extending the language: variant and inductive types

In this section, we show that the definition of forward AD and the semantics generalize if we extend the language of Section 3 with variants and inductive types. As an example of inductive types, we consider lists. This specific choice is only for expository purposes and the whole development works at the level of generality of arbitrary algebraic data types generated as initial algebras of (polynomial) type constructors formed by finite products and variants. These types are easily interpreted in the category of diffeological spaces in much the same way. The categorically minded reader may regard this as a consequence of $\mathbf{Diff}$ being a concrete Grothendieck quasitopos, e.g. [BH11], and hence is complete and cocomplete.

5.1. Language.

We additionally consider the following types and terms:

$\begin{array}[t]{l@{\quad\!\!}*3{l@{}}@{\,}l}{\tau},{\sigma},{\rho}&::=&\dots&% \mspace{-25.0mu}\qquad\text{types}\\ &\mathrel{\lvert}&\{\mathsf{\ell_{1}}\,{{\tau}_{1}}\mathrel{\big{\lvert}}% \ldots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{\tau}_{n}}\}&\qquad\text{% variant}\\ &\mathrel{\lvert}&\mathbf{list}({\tau})&\qquad\text{list}\\[6.0pt] {t},{s},{r}&::=&\dots&\mspace{-25.0mu}\qquad\text{terms}\\ &\mathrel{\lvert}&{\tau}.\ell\,{t}&\qquad\text{variant constructor}\\ &\mathrel{\lvert}&[\,]\ \mathrel{\lvert}\ {t}::{s}&\qquad\text{empty list and % cons}\\ &\mathrel{\lvert}&\mathbf{case}\,{t}\,\mathbf{of}\,\{\mathsf{\ell_{1}}\,{{x}_{% 1}}\to{{s}_{1}}\mathrel{\big{\lvert}}\cdots\mathrel{\big{\lvert}}\mathsf{\ell_% {n}}\,{{x}_{n}}\to{{s}_{n}}\}&\qquad\text{pattern matching: variants}\\ &\mathrel{\lvert}&\mathbf{fold}\,({x}_{1},{x}_{2}).{t}\,\mathbf{over}\,{s}\,% \mathbf{from}\,{r}&\qquad\text{list fold}\\ \end{array}$

We extend the type system according to the rules of Fig. 4.

$\begin{array}[]{@{}c@{}}\inferrule{\Gamma\vdash{t}:{\tau}_{i}}{\Gamma\vdash{% \tau}.\ell_{i}\,{t}:{\tau}}((\mathsf{\ell_{i}}\,{\tau}_{i})\in{\tau})\quad% \inferrule{~{}}{\Gamma\vdash[\,]:\mathbf{list}({\tau})}\quad\inferrule{\Gamma% \vdash{t}:{\tau}\\ \Gamma\vdash{s}:\mathbf{list}({\tau})}{\Gamma\vdash{t}::{s}:\mathbf{list}({% \tau})}\\ \\ \inferrule{\Gamma\vdash{t}:\{\mathsf{\ell_{1}}\,{{\tau}_{1}}\mathrel{\big{% \lvert}}\ldots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{\tau}_{n}}\}\\ \text{for each $1\leq i\leq n$: }\Gamma,{x}_{i}:{\tau}_{i}\vdash{s}_{i}:{\tau}% }{\Gamma\vdash\mathbf{case}\,{t}\,\mathbf{of}\,\{\begin{array}[t]{@{}l@{\,}l@{% }l@{}}\mathsf{\ell_{1}}\,{{x}_{1}}\to{{s}_{1}}\mathrel{\big{\lvert}}\cdots% \mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{x}_{n}}&\to{{s}_{n}}\}:{\tau}\end{% array}}\\ \\ \inferrule{\Gamma\vdash{s}:\mathbf{list}({\tau})\\ \Gamma\vdash{r}:{\sigma}\\ \Gamma,{x}_{1}:{\tau},{x}_{2}:{\sigma}\vdash{t}:{\sigma}}{\Gamma\vdash\mathbf{% fold}\,({x}_{1},{x}_{2}).{t}\,\mathbf{over}\,{s}\,\mathbf{from}\,{r}:{\sigma}}% \end{array}$

Figure 4. Additional typing rules for the extended language.

We can then extend $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ (again, writing it as $\overrightarrow{\mathcal{D}}$ , for legibility) to our new types and terms by

	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\{\mathsf{\ell_{1}% }\,{{\tau}_{1}}\mathrel{\big{\lvert}}\ldots\mathrel{\big{\lvert}}\mathsf{\ell_% {n}}\,{{\tau}_{n}}\})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\{\mathsf{\ell_% {1}}\,{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}_{1})}\mathrel{% \big{\lvert}}\ldots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({\tau}_{n})}\}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\mathbf{list}({% \tau}))\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{list}(\scalebox{0.8}{% $\overrightarrow{\mathcal{D}}$}({\tau}))$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({\tau}.\ell\,{t})% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}\scalebox{0.8}{$\overrightarrow{% \mathcal{D}}$}({\tau}).\ell\,{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({% t})}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}([\,])\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}[\,]$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({t}::{s})\stackrel% {{\scriptstyle\mathrm{def}}}{{=}}\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$% }({t})::\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}({s})$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\mathbf{case}\,{t}% \,\mathbf{of}\,\{\mathsf{\ell_{1}}\,{{x}_{1}}\to{{s}_{1}}\mathrel{\big{\lvert}% }\cdots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{x}_{n}}\to{{s}_{n}}\})% \stackrel{{\scriptstyle\mathrm{def}}}{{=}}$
	$\displaystyle\,\quad\mathbf{case}\,\scalebox{0.8}{$\overrightarrow{\mathcal{D}% }$}({t})\,\mathbf{of}\,\{\mathsf{\ell_{1}}\,{{x}_{1}}\to{\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({s}_{1})}\mathrel{\big{\lvert}}\cdots\mathrel{% \big{\lvert}}\mathsf{\ell_{n}}\,{{x}_{n}}\to{\scalebox{0.8}{$\overrightarrow{% \mathcal{D}}$}({s}_{n})}\}$
	$\displaystyle\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}(\mathbf{fold}\,({x% }_{1},{x}_{2}).{t}\,\mathbf{over}\,{s}\,\mathbf{from}\,{r})\stackrel{{% \scriptstyle\mathrm{def}}}{{=}}\mathbf{fold}\,({x}_{1},{x}_{2}).\scalebox{0.8}% {$\overrightarrow{\mathcal{D}}$}({t})\,\mathbf{over}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({s})\,\mathbf{from}\,\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}({r})$

To demonstrate the practical use of expressive type systems for differential programming, we consider the following two examples. {exa}[Lists of inputs for neural nets] Usually, we run a neural network on a large data set, the size of which might be determined at runtime. To evaluate a neural network on multiple inputs, in practice, one often sums the outcomes. This can be coded in our extended language as follows. Suppose that we have a network $f:\boldsymbol{(}\mathbf{real}^{n}\boldsymbol{\mathop{*}}P\boldsymbol{)}\to% \mathbf{real}$ that operates on single input vectors. We can construct one that operates on lists of inputs as follows:

g\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda\langle l,w\rangle.{\mathbf{% fold}\,({x}_{1},{x}_{2}).f\langle{x}_{1},w\rangle+{x}_{2}\,\mathbf{over}\,l\,% \mathbf{from}\,\underline{0}}:\boldsymbol{(}\mathbf{list}(\mathbf{real}^{n})% \boldsymbol{\mathop{*}}P\boldsymbol{)}\to\mathbf{real}

{exa}

[Missing data] In practically every application of statistics and machine learning, we face the problem of missing data: for some observations, only partial information is available.

In an expressive typed programming language like we consider, we can model missing data conveniently using the data type $\mathbf{maybe}({\tau})=\{\mathsf{\mathsf{Nothing}}\,{\boldsymbol{(}\,% \boldsymbol{)}}\mathrel{\big{\lvert}}\mathsf{\mathsf{Just}}\,{{\tau}}\}$ . In the context of a neural network, one might use it as follows. First, define some helper functions

		$\displaystyle\mathrm{fromMaybe}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \lambda{x}.{\lambda m.{\mathbf{case}\,m\,\mathbf{of}\,\{\mathsf{\mathsf{% Nothing}}\,{\_}\to{{x}}\mathrel{\big{\lvert}}\mathsf{\mathsf{Just}}\,{{x}^{% \prime}}\to{{x}^{\prime}}\}}}$
		$\displaystyle\mathrm{fromMaybe}^{n}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}% \lambda\langle{x}_{1},{.}{.}{.},{x}_{n}\rangle.{\lambda\langle m_{1},{.}{.}{.}% ,m_{n}\rangle.{\langle\mathrm{fromMaybe}\,{x}_{1}\,m_{1},{.}{.}{.},\mathrm{% fromMaybe}\,{x}_{n}\,m_{n}\rangle}}$
		$\displaystyle\qquad\qquad:(\mathbf{maybe}(\mathbf{real}))^{n}\to\mathbf{real}^% {n}\to\mathbf{real}^{n}$
		$\displaystyle\mathrm{map}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\lambda f.{% \lambda l.{\mathbf{fold}\,({x}_{1},{x}_{2}).f\,{x}_{1}::{x}_{2}\,\mathbf{over}% \,l\,\mathbf{from}\,[\,]}}:({\tau}\to{\sigma})\to\mathbf{list}({\tau})\to% \mathbf{list}({\sigma})$

Given a neural network $f:\boldsymbol{(}\mathbf{list}(\mathbf{real}^{k})\boldsymbol{\mathop{*}}P% \boldsymbol{)}\to\mathbf{real}$ , we can build a new one that operates on on a data set for which some covariates (features) are missing, by passing in default values to replace the missing covariates:

\lambda\langle l,\langle m,w\rangle\rangle.f\langle\mathrm{map}\,(\mathrm{% fromMaybe}^{k}\,m)\,l,w\rangle:\boldsymbol{(}\mathbf{list}((\mathbf{maybe}(% \mathbf{real}))^{k})\boldsymbol{\mathop{*}}\boldsymbol{(}\mathbf{real}^{k}% \boldsymbol{\mathop{*}}P\boldsymbol{)}\boldsymbol{)}\to\mathbf{real}

Then, given a data set $l$ with missing covariates, we can perform automatic differentiation on this network to optimize, simultaneously, the ordinary network parameters $w$ and the default values for missing covariates $m$ .

5.2. Semantics.

In Section 4 we gave a denotational semantics for the simple language in diffeological spaces. This extends to the language in this section, as follows. As before, each type ${\tau}$ is interpreted as a diffeological space, which is a set equipped with a family of plots:

•

A variant type $\{\mathsf{\ell_{1}}\,{{\tau}_{1}}\mathrel{\big{\lvert}}\ldots\mathrel{\big{% \lvert}}\mathsf{\ell_{n}}\,{{\tau}_{n}}\}$ is inductively interpreted as the disjoint union of the semantic spaces, $\textstyle\llbracket\{\mathsf{\ell_{1}}\,{{\tau}_{1}}\mathrel{\big{\lvert}}% \dots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{\tau}_{n}}\}\rrbracket\ \ % \stackrel{{\scriptstyle\mathrm{def}}}{{=}}\ \ \biguplus_{i=1}^{n}\llbracket{% \tau}_{i}\rrbracket$ , with $U$ -plots

\textstyle\mathcal{P}_{\llbracket\{\mathsf{\ell_{1}}\,{{\tau}_{1}}\mathrel{% \big{\lvert}}\ldots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{\tau}_{n}}\}% \rrbracket}^{U}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\{\left.\left[U_% {j}\xrightarrow{f_{j}}\llbracket{\tau}_{j}\rrbracket\to\biguplus_{i=1}^{n}% \llbracket{\tau}_{i}\rrbracket\right]_{j=1}^{n}~{}\right|~{}U=\biguplus_{j=1}^% {n}U_{j},\;f_{j}\in\mathcal{P}_{\llbracket{\tau}_{j}\rrbracket}^{U_{j}}\right\}.

•

A list type $\mathbf{list}({\tau})$ is interpreted as the union of the sets of length $i$ tuples for all natural numbers $i$ , $\llbracket\mathbf{list}({\tau})\rrbracket\ \ \stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\ \ \biguplus_{i=0}^{\infty}\llbracket{\tau}\rrbracket^{i}$ with $U$ -plots

\textstyle\mathcal{P}_{\llbracket\mathbf{list}({\tau})\rrbracket}^{U}\stackrel% {{\scriptstyle\mathrm{def}}}{{=}}\left\{\left.\left[U_{j}\xrightarrow{f_{j}}% \llbracket{\tau}\rrbracket^{j}\to\biguplus_{i=0}^{\infty}\llbracket{\tau}% \rrbracket^{i}\right]_{j=0}^{\infty}~{}\right|~{}U=\biguplus_{j=0}^{\infty}U_{% j},\;f_{j}\in\mathcal{P}_{\llbracket{\tau}\rrbracket^{j}}^{U_{j}}\right\}

The constructors and destructors for variants and lists are interpreted as in the usual set theoretic semantics.

It is routine to show inductively that these interpretations are smooth. Thus every term $\Gamma\vdash{t}:{\tau}$ in the extended language is interpreted as a smooth function $\llbracket{t}\rrbracket:\llbracket\Gamma\rrbracket\to\llbracket{\tau}\rrbracket$ between diffeological spaces. List objects as initial algebras are computed as usual in a cocomplete category (e.g. [JR11]). More generally, the interpretation for algebraic data types follows exactly the usual categorical semantics of variant types and inductive types (e.g. [Pit95]).

6. Categorical analysis of (higher order) forward AD and its correctness

This section has three parts. First, we give a categorical account of the functoriality of AD (Ex. 6.1). Then we introduce our gluing construction, and relate it to the correctness of AD (dgm. (4)). Finally, we state and prove a correctness theorem for all first order types by considering a category of manifolds (Th. 8).

6.1. Syntactic categories.

The key contribution of this subsection is that the AD macro translation (Section 3.2) has a canonical status as a unique functor between categories with structure. To this end, we build a syntactic category $\mathbf{Syn}$ from our language, which has the property of a free category with certain structure. This means that for any category $\mathcal{C}$ with this structure, there is a unique structure-preserving functor $\mathbf{Syn}\to\mathcal{C}$ , which is an interpretation of our language in that category. Generally speaking, this is the categorical view of denotational semantics (e.g. [Pit95]). But in this particular setting, the category $\mathbf{Syn}$ itself admits alternative forms of this structure, given by the dual numbers interpretation, the triple numbers interpretation etc. of Section 2. This gives canonical functors $\mathbf{Syn}\to\mathbf{Syn}$ translating the language into itself, which are the AD macro translations (Section 3.2). A key point is that $\mathbf{Syn}$ is almost entirely determined by universal properties (for example, cartesian closure for the function space); the only freedom is in the choice of interpretation of

(1)

the real numbers $\mathbf{real}$ , which can be taken as the plain type $\mathbf{real}$ , or as the dual numbers interpretation $\mathbf{real}\ast\mathbf{real}$ etc.;
(2)

the primitive operations $\mathsf{op}$ , which can be taken as the operation $\mathsf{op}$ itself, or as the derivative of the operation etc..

$\displaystyle\mathbf{case}\,\langle{t}_{1},\ldots,{t}_{n}\rangle\,\mathbf{of}% \,\langle{x}_{1},\ldots,{x}_{n}\rangle\to{s}={s}{}[^{{t}_{1}}\!/\!_{{x}_{1}},% \ldots,^{{t}_{n}}\!/\!_{{x}_{n}}]$ $\displaystyle{s}{}[^{{t}}\!/\!_{{y}}]\stackrel{{\scriptstyle\#{x}_{1},\ldots,{% x}_{n}}}{{=}}\mathbf{case}\,{t}\,\mathbf{of}\,\langle{x}_{1},\ldots,{x}_{n}% \rangle\to{s}{}[^{\langle{x}_{1},\ldots,{x}_{n}\rangle}\!/\!_{{y}}]$ $\displaystyle\mathbf{case}\,\mathsf{\ell_{i}}\,{{t}}\,\mathbf{of}\,\{\mathsf{% \ell_{1}}\,{{x}_{1}}\to{{s}_{1}}\mathrel{\big{\lvert}}\cdots\mathrel{\big{% \lvert}}\mathsf{\ell_{n}}\,{{x}_{n}}\to{{s}_{n}}\}={s}_{i}{}[^{{t}}\!/\!_{{x}_% {i}}]$ $\displaystyle{s}{}[^{{t}}\!/\!_{{y}}]\stackrel{{\scriptstyle\#{x}_{1},\ldots,{% x}_{n}}}{{=}}\mathbf{case}\,{t}\,\mathbf{of}\,\{\mathsf{\ell_{1}}\,{{x}_{1}}% \to{{s}{}[^{\mathsf{\ell_{1}}\,{{x}_{1}}}\!/\!_{{y}}]}\mathrel{\big{\lvert}}% \cdots\mathrel{\big{\lvert}}\mathsf{\ell_{n}}\,{{x}_{n}}\to{{s}{}[^{\mathsf{% \ell_{n}}\,{{x}_{n}}}\!/\!_{{y}}]}\}$ $\displaystyle\mathbf{fold}\,({x}_{1},{x}_{2}).{t}\,\mathbf{over}\,[\,]\,% \mathbf{from}\,{r}={r}$ $\displaystyle\mathbf{fold}\,({x}_{1},{x}_{2}).{t}\,\mathbf{over}\,{s}_{1}::{s}% _{2}\,\mathbf{from}\,{r}={t}{}[^{{s}_{1}}\!/\!_{{x}_{1}},^{\mathbf{fold}\,({x}% _{1},{x}_{2}).{t}\,\mathbf{over}\,{s}_{2}\,\mathbf{from}\,{r}}\!/\!_{{x}_{2}}]$ $\displaystyle u={s}{}[^{[\,]}\!/\!_{{y}}],{r}{}[^{{s}}\!/\!_{{x}_{2}}]={s}{}[^% {{x}_{1}::{y}}\!/\!_{{y}}]\Rightarrow{s}{}[^{{t}}\!/\!_{{y}}]\stackrel{{% \scriptstyle\#{x}_{1},{x}_{2}}}{{=}}\mathbf{fold}\,({x}_{1},{x}_{2}).{r}\,% \mathbf{over}\,{t}\,\mathbf{from}\,u$ $\displaystyle(\lambda{x}.{{t}})\,{s}={t}{}[^{{s}}\!/\!_{{x}}]$ $\displaystyle{t}\stackrel{{\scriptstyle\#{x}}}{{=}}\lambda{x}.{{t}\,{x}}$ We write $\stackrel{{\scriptstyle\#{x}_{1},\ldots,{x}_{n}}}{{=}}$ to indicate that the variables are free in the left hand side

Figure 5. Standard

\beta\eta

-laws (e.g. [Pit95]) for products, functions, variants and lists.

In more detail, our language induces a syntactic category as follows. {defi} Let $\mathbf{Syn}$ be the category whose objects are types, and where a morphism ${\tau}\to{\sigma}$ is a term in context ${x}:{\tau}\vdash{t}:{\sigma}$ modulo the $\beta\eta$ -laws (Fig. 5). Composition is by substitution. For simplicity, we do not impose identities involving the primitive operations, such as the arithmetic identity $x+y=y+x$ in $\mathbf{Syn}$ . As is standard, this category has the following universal property.

Lemma 4 (e.g. [Pit95]).

For every bicartesian closed category $\mathcal{C}$ with list objects, and every choice of an object $F(\mathbf{real})\in\mathcal{C}$ and morphisms $F(\mathsf{op})\in\mathcal{C}(F(\mathbf{real})^{n},F(\mathbf{real}))$ for all $\mathsf{op}\in\mathsf{Op}_{n}$ and $n\in\mathbb{N}$ , in $\mathcal{C}$ , there is a unique functor $F:{\mathbf{Syn}\to\mathcal{C}}$ respecting the interpretation and preserving the bicartesian closed structure as well as list objects.

Proof 6.1 (Proof notes).

The functor $F:\mathbf{Syn}\to\mathcal{C}$ is a canonical denotational semantics for the language, interpreting types as objects of $\mathcal{C}$ and terms as morphisms. For instance, $F({{\tau}\to{\sigma}})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(F{\tau}\to F{% {\sigma}})$ , the function space in the category $\mathcal{C}$ , and $F{({t}\,{s})}$ is the composite $(F{{t}},F{{s}});\mathit{eval}$ .

When $\mathcal{C}=\mathbf{Diff}$ , the denotational semantics of the language in diffeological spaces (Section 4,5.2) can be understood as the unique structure preserving functor $\llbracket-\rrbracket:\mathbf{Syn}\to\mathbf{Diff}$ satisfying $\llbracket\mathbf{real}\rrbracket=\mathbb{R}$ , $\llbracket\varsigma\rrbracket=\varsigma$ and so on.

{exa}

[Canonical definition of forward AD] The forward AD macro $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ (Section 3,5.1) arises as a canonical bicartesian closed functor on $\mathbf{Syn}$ that preserves list objects. Consider the unique bicartesian closed functor $F:\mathbf{Syn}\to\mathbf{Syn}$ that preserves list objects such that $F(\mathbf{real})=\mathbf{real}^{\binom{R+k}{k}}$ and

F(\mathsf{op})={z}\!:\!\boldsymbol{(}F(\mathbf{real})\boldsymbol{\mathop{*}}..% \boldsymbol{\mathop{*}}F(\mathbf{real})\boldsymbol{)}\vdash\mathbf{case}\,{z}% \,\mathbf{of}\,\langle{x}_{1},...,{x}_{n}\rangle\to\scalebox{0.8}{$% \overrightarrow{\mathcal{D}}$}_{\!(k,R)}(\mathsf{op}({x}_{1},\ldots,{x}_{n})):% F(\mathbf{real}).

Then for any type ${\tau}$ , $F({\tau})=\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({\tau})$ , and for any term $x:{\tau}\vdash{t}:{\sigma}$ , $F({t})=\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})$ as morphisms $F({\tau})\to F({\sigma})$ in the syntactic category.

This observation is a categorical counterpart to Lemma 1.

6.2. Categorical gluing and logical relations.

Gluing is a method for building new categorical models which has been used for many purposes, including logical relations and realizability [MS92]. Our logical relations argument in the proof of Theorem 3 can be understood in this setting. (In fact we originally found the proof of Theorem 3 in this way.) In this subsection, for the categorically minded, we explain this, and in doing so we quickly recover a correctness result for the more general language in Section 5 and for arbitrary first order types.

The general, established idea of categorical logical relations starts from the observation that that logical relations are defined by induction on the structure of types. Types have universal properties in a categorical semantics (e.g. cartesian closure for the function space), and so we can organize the logical relations argument by defining some category $\mathcal{C}$ of relations and observing that it has the requisite categorical structure. The interpretation of types as relations can then be understood as coming from a unique structure preserving map $\mathbf{Syn}\to\mathcal{C}$ . In this paper, our logical relations are not quite as simple as a binary relation on sets; rather they are relations between plots. Nonetheless, this still forms a category with the appropriate structure, which follows because it can still be regarded as arising from a gluing construction, as we now explain.

We define a category $\mathbf{Gl}_{k}$ whose objects are triples $(X,X^{\prime},S)$ where $X$ and $X^{\prime}$ are diffeological spaces and $S\subseteq\mathcal{P}_{X}^{\mathbb{R}^{k}}\times\mathcal{P}_{X^{\prime}}^{% \mathbb{R}^{k}}$ is a relation between their $k$ -dimensional plots. A morphism $(X,X^{\prime},S)\to(Y,Y^{\prime},T)$ is a pair of smooth functions $f\colon X\to Y$ , $f^{\prime}\colon X^{\prime}\to Y^{\prime}$ , such that if $(g,g^{\prime})\in S$ then $(g;f,g^{\prime};f^{\prime})\in T$ . The idea is that this is a semantic domain where we can simultaneously interpret the language and its automatic derivatives.

Proposition 5.

The category $\mathbf{Gl}_{k}$ is bicartesian closed, has list objects, and the projection functor $\mathrm{proj}:\mathbf{Gl}_{k}\to\mathbf{Diff}\times\mathbf{Diff}$ preserves this structure.

Proof 6.2 (Proof notes).

The category $\mathbf{Gl}_{k}$ is a full subcategory of the comma category
${\rm id}_{\mathbf{Set}}\downarrow\mathbf{Diff}(\mathbb{R}^{k},-)\times\mathbf{% Diff}(\mathbb{R}^{k},-)$ . The result thus follows by the general theory of categorical gluing (e.g. [JLS07, Lemma 15]).

We give a semantics $\llparenthesis-\rrparenthesis=(\llparenthesis-\rrparenthesis_{0},% \llparenthesis-\rrparenthesis_{1},S_{-})$ for the language in $\mathbf{Gl}_{k}$ , interpreting types ${\tau}$ as objects $(\llparenthesis{\tau}\rrparenthesis_{0},\llparenthesis{\tau}\rrparenthesis_{1}% ,S_{{\tau}})$ , and terms as morphisms. We let $\llparenthesis\mathbf{real}\rrparenthesis_{0}\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\mathbb{R}$ and $\llparenthesis\mathbf{real}\rrparenthesis_{1}\stackrel{{\scriptstyle\mathrm{% def}}}{{=}}\mathbb{R}^{\binom{R+k}{k}}$ , with the relation

S_{\mathbf{real}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\{(f,{\left(% \frac{\partial^{\alpha_{1}+\ldots+\alpha_{k}}f(x)}{\partial x_{1}^{\alpha_{1}}% \cdots\partial x_{k}^{\alpha_{k}}}\right)}_{(\alpha_{1},...,\alpha_{k})=(0,...% ,0)}^{(R,0,...,0)})~{}|~{}f:\mathbb{R}^{k}\to\mathbb{R}\text{ smooth}\right\}.

We interpret the operations $\mathsf{op}$ according to $\llbracket\mathsf{op}\rrbracket$ in $\llparenthesis-\rrparenthesis_{0}$ , but according to the $(k,R)$ -Taylor representation of $\llbracket\mathsf{op}\rrbracket$ in $\llparenthesis-\rrparenthesis_{1}$ . For instance, when $k=2$ and $r=2$ , $\llparenthesis*\rrparenthesis_{1}:\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb% {R}^{2}$ is

	$\displaystyle\llparenthesis*\rrparenthesis_{1}$	$\displaystyle((x_{00},x_{01},x_{02},x_{10},x_{11},x_{20}),(y_{00},y_{01},y_{02% },y_{10},y_{11},y_{20}))\stackrel{{\scriptstyle\mathrm{def}}}{{=}}$
		$\displaystyle(x_{00}y_{00},$
		$\displaystyle\;x_{00}y_{01}+x_{01}y_{00},$
		$\displaystyle\;x_{02}y_{00}+2x_{01}y_{01}+x_{00}y_{02},$
		$\displaystyle\;x_{00}y_{10}+x_{10}y_{00},$
		$\displaystyle\;x_{11}y_{00}+x_{01}y_{10}+x_{10}y_{01}+x_{00}y_{11},$
		$\displaystyle\;x_{20}y_{00}+2x_{10}y_{10}+x_{00}y_{20})\text{.}$

At this point one checks that these interpretations are indeed morphisms in $\mathbf{Gl}_{k}$ . This is equivalent to the statement that $\llparenthesis\mathsf{op}\rrparenthesis_{1}$ is the $(k,R)$ -Taylor representation of $\llbracket\mathsf{op}\rrbracket$ (3). The remaining constructions of the language are interpreted using the categorical structure of $\mathbf{Gl}_{k}$ , following Lemma 4.

Notice that the diagram below commutes. One can check this by hand or note that it follows from the initiality of $\mathbf{Syn}$ (Lemma 4): all the functors preserve all the structure.

(4)

We thus arrive at a restatement of the correctness theorem (Th. 3), which holds even for the extended language with variants and lists, because for any $x_{1}:\mathbf{real},{.}{.}{.},x_{n}:\mathbf{real}\vdash{t}:\mathbf{real}$ , the interpretations $(\llbracket{t}\rrbracket,\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}% }$}_{(k,R)}({t})\rrbracket)$ are in the image of the projection $\mathbf{Gl}_{k}\to\mathbf{Diff}\times\mathbf{Diff}$ , and hence $\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({t})\rrbracket$ is a $(k,R)$ -Taylor representation of $\llbracket{t}\rrbracket$ .

6.3. Correctness at all first order types, via manifolds.

We now generalize Theorem 3 to hold at all first order types, not just the reals.

So far, we have shown that our macro translation (Section 3.2) gives correct derivatives to functions of the real numbers, even if other types are involved in the definitions of the functions (Theorem 3 and Section 6.2). We can state this formally because functions of the real numbers have well understood derivatives (Section 2). There are no established mathematical notions of derivatives at higher types, and so we cannot even begin to argue that our syntactic derivatives of functions $(\mathbf{real}\to\mathbf{real})\to(\mathbf{real}\to\mathbf{real})$ match with some existing mathematical notion (see also Section 7).

However, for functions of first order type, like $\mathbf{list}(\mathbf{real})\to\mathbf{list}(\mathbf{real})$ , there are established mathematical notions of derivative, because we can understand $\mathbf{list}(\mathbf{real})$ as the manifold of all tuples of reals, and then appeal to the well-known theory of manifolds and jet bundles. We do this now, to achieve a correctness theorem for all first order types (Theorem 8). The key high level points are that

•

manifolds support a notion of differentiation, and an interpretation of all first order types, but not an interpretation of higher types;
•

diffeological spaces support all types, including higher types, but not an established notion of differentiation in general;
•

manifolds and smooth maps embed full and faithfully in diffeological spaces, preserving the interpretation of first order types, so we can use the two notions together.

We now explain this development in more detail.

For our purposes, a smooth manifold $M$ is a second-countable Hausdorff topological space together with a smooth atlas. In more detail, a topological space $X$ is second-countable when there exists a collection $U:=\{U_{i}\}_{i\in\mathbb{N}}$ of open subsets of $X$ such that any open subset of $X$ can be written as a union of elements of $U$ . A topological space $X$ is Hausdorff if for every distinct points $x$ and $y$ , there exists disjoint open subsets $U,V$ of $X$ such that $x\in U,y\in V$ . A smooth atlas of a topological space $X$ is an open cover $\mathcal{U}$ together with homeomorphisms $\left(\phi_{U}:U\to\mathbb{R}^{n(U)}\right)_{U\in\mathcal{U}}$ (called charts, or local coordinates) such that $\phi_{U}^{-1};\phi_{V}$ is smooth on its domain of definition for all $U,V\in\mathcal{U}$ . A function $f:M\to N$ between manifolds is smooth if $\phi^{-1}_{U};f;\psi_{V}$ is smooth for all charts $\phi_{U}$ and $\psi_{V}$ of $M$ and $N$ , respectively. Let us write $\mathbf{Man}$ for this category. This definition of manifolds is a slight generalisation of the more usual one from differential geometry because different charts in an atlas may have different finite dimensions $n(U)$ . Thus we consider manifolds with dimensions that are potentially unbounded, albeit locally finite.

Each open subset of $\mathbb{R}^{n}$ can be regarded as a manifold. This lets us regard the category of manifolds $\mathbf{Man}$ as a full subcategory of the category of diffeological spaces. We consider a manifold $(X,\{\phi_{U}\}_{U})$ as a diffeological space with the same carrier set $X$ and where the plots $\mathcal{P}_{X}^{U}$ , called the manifold diffeology, are the smooth functions in $\mathbf{Man}(U,X)$ . A function $X\to Y$ is smooth in the sense of manifolds if and only if it is smooth in the sense of diffeological spaces [IZ13]. For the categorically minded reader, this means that we have a full embedding of $\mathbf{Man}$ into $\mathbf{Diff}$ . Moreover, the natural interpretation of the first order fragment of our language in $\mathbf{Man}$ coincides with that in $\mathbf{Diff}$ . That is, the embedding of $\mathbf{Man}$ into $\mathbf{Diff}$ preserves finite products and countable coproducts (hence initial algebras of polynomial endofunctors).

Proposition 6.

Suppose that a type ${\tau}$ is first order, i.e. it is just built from reals, products, variants, and lists (or, again, arbitrary inductive types), and not function types. Then the diffeological space $\llbracket{\tau}\rrbracket$ is a manifold.

Proof 6.3 (Proof notes).

This is proved by induction on the structure of types. In fact, one may show that every such $\llbracket{\tau}\rrbracket$ is isomorphic to a manifold of the form $\biguplus_{i=1}^{n}\mathbb{R}^{d_{i}}$ where the bound $n$ is either finite or $\infty$ , but this isomorphism is typically not an identity function.

We recall how the Taylor representation of any morphism $f:M\to N$ of manifolds is given by its action on jets [KSM99, Chapter IV]. For each point $x$ in a manifold $M$ , define the $(k,R)$ -jet space $\mathcal{J}^{(k,R)}_{x}M$ to be the set $\{\gamma\in\mathbf{Man}(\mathbb{R}^{k},M)\mid\gamma(0)=x\}/\sim$ of equivalence classes $[\gamma]$ of $k$ -dimensional plots $\gamma$ in $M$ based at $x$ , where we identify $\gamma_{1}\sim\gamma_{2}$ iff all partial derivatives of order $\leq R$ coincide in the sense that

\frac{\partial^{\alpha_{1}+...+\alpha_{k}}(\gamma_{1};f)(x)}{\partial x_{1}^{% \alpha_{1}}\cdots\partial x_{k}^{\alpha_{k}}}(0)=\frac{\partial^{\alpha_{1}+..% .+\alpha_{k}}(\gamma_{2};f)(x)}{\partial x_{1}^{\alpha_{1}}\cdots\partial x_{k% }^{\alpha_{k}}}(0)

for all smooth $f:M\to\mathbb{R}$ and all multi-indices $(\alpha_{1},...,\alpha_{k})=(0,...,0),...,(R,0,...,0)$ . In the case of $(k,R)=(1,1)$ , a $(k,R)$ -jet space is better known as a tangent space. The $(k,R)$ -jet bundle (a.k.a. tangent bundle, in case $(k,R)=(1,1)$ ) of $M$ is the set $\mathcal{J}^{(k,R)}(M)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\biguplus_{x% \in M}\mathcal{J}^{(k,R)}_{x}(M)$ . The charts of $M$ equip $\mathcal{J}^{(k,R)}(M)$ with a canonical manifold structure. The (manifold) diffeology of these jet bundles can be concisely summarized by the plots $\mathcal{P}_{\mathcal{J}^{(k,R)}(M)}^{U}=\left\{f:U\to|\mathcal{J}^{(k,R)}(M)|% \mid\exists g\in\mathcal{P}_{M}^{U\times\mathbb{R}^{k}}.\forall u\in U.(g(u,0)% ,[v\mapsto g(u,v)])=f(u)\right\}$ .
Then $\mathcal{J}^{(k,R)}$ acts on smooth maps $f:M\to N$ to give $\mathcal{J}^{(k,R)}(f):\mathcal{J}^{(k,R)}(M)\to\mathcal{J}^{(k,R)}(N)$ is defined as $\mathcal{J}^{(k,R)}(f)(x,[\gamma])\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(f% (x),[\gamma;f])$ . In local coordinates, this action $\mathcal{J}^{(k,R)}(f)$ is seen to coincide precisely with the $(k,R)$ -Taylor representation of $f$ given by the Faà di Bruno formula [Mer04]. All told, the $(k,R)$ -jet bundle is a functor $\mathcal{J}^{(k,R)}:\mathbf{Man}\to\mathbf{Man}$ [KSM99].

We can understand the jet bundle of a composite space in terms of that of its parts.

Lemma 7.

There are canonical isomorphisms $\mathcal{J}^{(k,R)}(\biguplus_{i=1}^{\infty}M_{i})\cong\biguplus_{i=1}^{\infty% }\mathcal{J}^{(k,R)}(M_{i})$ and $\mathcal{J}^{(k,R)}(M_{1}\times\ldots\times M_{n})\cong\mathcal{J}^{(k,R)}(M_{% 1})\times\ldots\times\mathcal{J}^{(k,R)}(M_{n})$ .

Proof 6.4 (Proof notes).

For disjoint unions, notice that that smooth morphisms from $\mathbb{R}^{k}$ into a disjoint union of manifolds always factor over a single inclusion, because $\mathbb{R}^{k}$ is connected. For products, it is well-known that partial derivatives of a morphism $(f_{1},...,f_{n})$ are calculated component-wise [Lee13, ex. 3-2].

We define a canonical isomorphism $\phi_{{\tau}}^{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}\mathcal{J}}:% \llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}({\tau})% \rrbracket\to\mathcal{J}^{(k,R)}(\llbracket{\tau}\rrbracket)$ for every type ${\tau}$ , by induction on the structure of types. We let $\phi_{\mathbf{real}}^{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}\mathcal{J% }}:\llbracket\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}(\mathbf{% real})\rrbracket\to\mathcal{J}^{(k,R)}(\llbracket\mathbf{real}\rrbracket)$ be given by $\phi_{\mathbf{real}}^{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}\mathcal{J% }}(x,x^{\prime})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(x,[t\mapsto x+x^{% \prime}t])$ . For the other types, we use Lemma 7. We can now phrase correctness at all first order types.

Theorem 8 (Semantic correctness of ${\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}}$ (full)).

For any ground ${\tau}$ , any first order context $\Gamma$ and any term $\Gamma\vdash{t}:{\tau}$ , the syntactic translation $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}$ coincides with the $(k,R)$ -jet bundle functor, modulo these canonical isomorphisms:

Proof 6.5 (Proof notes).

For any $k$ -dimensional plot $\gamma\in\mathbf{Man}(\mathbb{R}^{k},M)$ , let $\bar{\gamma}\in\mathbf{Man}(\mathbb{R}^{k},\mathcal{J}^{(k,R)}(M))$ be the $(k,R)$ -jet curve, given by $\bar{\gamma}(x)=(\gamma(x),[t\mapsto\gamma(x+t)])$ . First, we note that a smooth map $h:\mathcal{J}^{(k,R)}(M)\to\mathcal{J}^{(k,R)}(N)$ is of the form $\mathcal{J}^{(k,R)}(g)$ for some $g:M\to N$ if for all smooth $\gamma:\mathbb{R}^{k}\to M$ we have $\bar{\gamma};h=\overline{(\gamma;g)}:\mathbb{R}^{k}\to\mathcal{J}^{(k,R)}(N)$ . This generalizes (3). Second, for any first order type ${\tau}$ , $S_{\llbracket{\tau}\rrbracket}=\{(f,\tilde{f})~{}|~{}\tilde{f};\phi_{\tau}^{% \scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}\mathcal{J}}=\bar{f}\}$ . This is shown by induction on the structure of types. We conclude the theorem from diagram (4), by putting these two observations together.

7. Discussion: What are derivatives of higher order functions?

In our gluing categories $\mathbf{Gl}_{k}$ of Section 6.2, we have avoided the question of what semantic derivatives should be associated with higher order functions. Our syntactic macro $\overrightarrow{\mathcal{D}}$ provides a specific derivative for every definable function, but in the model $\mathbf{Gl}_{k}$ there is only a relation between plots and their corresponding Taylor representations, and this relation is not necessarily single-valued. Our approach has been rather indifferent about what “the” correct derivative of a higher order function should be. Instead, all we have cared about is that we are using “a” derivative that is correct in the sense that it can never be used to produce incorrect derivatives for first order functions, where we do have an unambiguous notion of correct derivative.

7.1. Automatic derivatives of higher order functions may not be unique!

For a concrete example to show that derivatives of higher order functions might not be unique in our framework, let us consider the case $(k,R)=(1,1)$ and focus on first derivatives of the evaluation function

	$\displaystyle\mathrm{ev}:\$	$\displaystyle\mathbb{R}\to\llbracket(\mathbf{real}\to\mathbf{real})\to\mathbf{% real}\rrbracket=\mathbb{R}\to(\mathbb{R}\Rightarrow\mathbb{R})\Rightarrow% \mathbb{R};$
		$\displaystyle r\mapsto(f\mapsto f(r)).$

Our macro $\overrightarrow{\mathcal{D}}$ will return $\lambda a:\mathbb{R}.\lambda f:\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R}.f(a,1)$ . In this section we show that the lambda term $\lambda a:\mathbb{R}.\lambda f:\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R}.\mathrm{sort}f(a,1)$ is also a valid derivative of the evaluation map, where $\mathrm{sort}:(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}\times\mathbb{R}% )\Rightarrow(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}\times\mathbb{R})$ is defined by

\displaystyle\mathrm{sort}:=\lambda f.\lambda((r,\_).(\pi_{1}(f(r,0)),\pi_{2}% \scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}((-,0);f;\pi_{1})(r))).

This map is idempotent and it converts any map $\mathbb{R}\times\mathbb{R}\to\mathbb{R}\times\mathbb{R}$ into the dual-numbers representation of its first component. For example, $(\mathrm{sort}(\mathrm{swap}))$ is the constantly $(0,0)$ function, where we write

	$\displaystyle\mathrm{swap}:\$	$\displaystyle\mathbb{R}\times\mathbb{R}\to\mathbb{R}\times\mathbb{R}$
		$\displaystyle(r,r^{\prime})\mapsto(r^{\prime},r).$

According to our gluing semantics, a function $g:\llparenthesis{\tau}\rrparenthesis_{1}\to\llparenthesis{\sigma}% \rrparenthesis_{1}$ defines a correct $(k,R)$ -Taylor representation of a function $f:\llparenthesis{\tau}\rrparenthesis_{0}\to\llparenthesis{\sigma}% \rrparenthesis_{0}$ iff $(f,g)$ defines a morphism $\llparenthesis{\tau}\rrparenthesis\to\llparenthesis{\sigma}\rrparenthesis$ in $\mathbf{Gl}_{k}$ . In particular, there is no guarantee that every $f$ has a unique correct $(k,R)$ -Taylor representation $g$ . (Although such Taylor representations are, in fact, unique when ${\tau},{\sigma}$ are first order types.) The gluing relation $\llparenthesis(\mathbf{real}\to\mathbf{real})\to\mathbf{real}\rrparenthesis$ in $\mathbf{Gl}_{1}$ relates curves in $\gamma:\mathbb{R}\to(\mathbb{R}\Rightarrow\mathbb{R})\Rightarrow\mathbb{R}$ to “tangent curves” $\gamma^{\prime}:\mathbb{R}\to(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R})\Rightarrow\mathbb{R}\times\mathbb{R}$ . In this relation, the function $\mathrm{ev}$ is related to at least two different tangent curves.

Lemma 9.

We have a smooth map

	$\displaystyle\mathrm{sort}:\$	$\displaystyle(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}\times\mathbb{R})% \to(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}\times\mathbb{R})$
		$\displaystyle f\mapsto((r,\_)\mapsto(\pi_{1}(f(r,0)),\nabla((-,0);f;\pi_{1})(r% ))).$

Proof 7.1.

Let $f\in\mathcal{P}_{\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}\times\mathbb{% R}}^{U}$ and let $\gamma_{1},\gamma_{2}\in\mathcal{P}_{\mathbb{R}}^{U}$ . Then, also $u\mapsto(\gamma_{1}(u),0);f(u);\pi_{1}\in\mathcal{P}_{\mathbb{R}}^{U}=\mathbf{% Man}(U,\mathbb{R})$ by definition of the exponential in $\mathbf{Diff}$ . Therefore, we also have that $u\mapsto\nabla(u^{\prime}\mapsto(\gamma_{1}(u^{\prime}),0);f(u^{\prime});\pi_{% 1})\in\mathbf{Man}(U,\mathbb{R})=\mathcal{P}_{\mathbb{R}}^{U}$ , as we are working with infinitely differentiable smooth maps. Consequently,

u\mapsto(f;\mathrm{sort})(u)(\gamma_{1}(u),\gamma_{2}(u))=(\pi_{1}(f(u))(% \gamma_{1}(u),0),\nabla(u^{\prime}\mapsto(\gamma_{1}(u^{\prime}),0);f(u^{% \prime});\pi_{1}))\in\mathcal{P}_{\mathbb{R}\times\mathbb{R}}^{U},

by definition of the product in $\mathbf{Diff}$ . It follows that $(f;\mathrm{sort})\in\mathcal{P}_{\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{% R}\times\mathbb{R}}^{U}$ .

Proposition 10.

We have that both $(\mathrm{ev},\mathrm{ev}^{\prime}_{1})\in\llparenthesis(\mathbf{real}\to% \mathbf{real})\to\mathbf{real}\rrparenthesis$ and $(\mathrm{ev},\mathrm{ev}^{\prime}_{2})\in\llparenthesis(\mathbf{real}\to% \mathbf{real})\to\mathbf{real}\rrparenthesis$ for

	$\displaystyle\mathrm{ev}^{\prime}_{1}:$	$\displaystyle\mathbb{R}\to(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R})\Rightarrow\mathbb{R}\times\mathbb{R}$
		$\displaystyle a\mapsto(f\mapsto f(a,1))$
	$\displaystyle\mathrm{ev}^{\prime}_{2}:$	$\displaystyle\mathbb{R}\to(\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R})\Rightarrow\mathbb{R}\times\mathbb{R}$
		$\displaystyle a\mapsto(f\mapsto(\mathrm{sort}f)(a,1)).$

Proof 7.2.

By definition of $\llparenthesis-\rrparenthesis$ , we need to show that for any $(\gamma,\gamma^{\prime})\in\llparenthesis\mathbf{real}\to\mathbf{real}\rrparenthesis$ , we have that $(x\mapsto\mathrm{ev}(x)(\gamma(x)),x\mapsto\mathrm{ev}^{\prime}_{i}(x)(\gamma^% {\prime}(x)))\in\llparenthesis\mathbf{real}\rrparenthesis$ . This means that we need to show that for $i=1,2$

\displaystyle x\mapsto\mathrm{ev}^{\prime}_{i}(x)(\gamma^{\prime}(x))=(x% \mapsto\mathrm{ev}(x)(\gamma(x)),\nabla(x\mapsto\mathrm{ev}(x)(\gamma(x))))

Unrolling further, this means we need to show that for any $\gamma:\mathbb{R}\to\mathbb{R}\Rightarrow\mathbb{R}$ and $\gamma^{\prime}:\mathbb{R}\to\mathbb{R}\times\mathbb{R}\Rightarrow\mathbb{R}% \times\mathbb{R}$ such that for any $(\delta,\delta^{\prime})\in\llparenthesis\mathbf{real}\rrparenthesis$ (which means that $\delta:\mathbb{R}\to\mathbb{R}$ and $\delta^{\prime}=(\delta,\nabla\delta)$ ), we have that

\displaystyle\Big{(}r\mapsto\gamma(r)(\delta(r)),r\mapsto\gamma^{\prime}(r)(% \delta^{\prime}(r))\Big{)}\in\llparenthesis\mathbf{real}\rrparenthesis

The latter part finally means that we need to show that

\displaystyle r\mapsto\gamma^{\prime}(r)(\delta(r),\nabla\delta(r))=(r\mapsto% \gamma(r)(\delta(r)),\nabla(r\mapsto\gamma(r)(\delta(r))))

Now, focussing on $\mathrm{ev}^{\prime}_{1}$ : we need to show that

	$\displaystyle x\mapsto\mathrm{ev}^{\prime}_{1}(x)(\gamma^{\prime}(x))=(x% \mapsto\mathrm{ev}(x)(\gamma(x)),$
	$\displaystyle\nabla(x\mapsto\mathrm{ev}(x)(\gamma(x))))$

Inlining the definition of $\mathrm{ev}^{\prime}_{1}$ : we need to show that

\displaystyle x\mapsto\gamma^{\prime}(x)(x,1)=(x\mapsto\gamma(x)(x),\nabla(x% \mapsto\gamma(x)(x)))

This follows by assumption by choosing $\delta(r)=r$ , and hence $\delta^{\prime}(r)=(r,1)$ .

Focussing on $\mathrm{ev}^{\prime}_{2}$ : we need to show that

\displaystyle x\mapsto\mathrm{ev}^{\prime}_{2}(x)(\gamma^{\prime}(x))=(x% \mapsto\mathrm{ev}(x)(\gamma(x)),\nabla(x\mapsto\mathrm{ev}(x)(\gamma(x))))

Inlining $\mathrm{ev}^{\prime}_{2}$ ’s definition: we need to show that

	$\displaystyle\Big{(}x\mapsto((\pi_{1}(\gamma^{\prime}(x)(x,0)),x\mapsto\nabla(% (-,0);\gamma^{\prime}(x);\pi_{1})(x)))\Big{)}=$
	$\displaystyle\Big{(}x\mapsto((r,\_)\mapsto(\pi_{1}(\gamma^{\prime}(x)(r,0)),% \nabla((-,0);\gamma^{\prime}(x);\pi_{1})(r)))(x,1)\Big{)}$

is equal to

	$\displaystyle x\mapsto(\mathrm{sort}~{}\gamma^{\prime}(x))(x,1)=$
	$\displaystyle\Big{(}x\mapsto\gamma(x)(x),\nabla(x\mapsto\gamma(x)(x))\Big{)}$

That is, we need to show that $\pi_{1}(\gamma^{\prime}(x)(x,0))=\gamma(x)(x)$ for all $x\in\mathbb{R}$ , which holds by the assumption that $(\gamma,\gamma^{\prime})\in\llparenthesis\mathbf{real}\to\mathbf{real}\rrparenthesis$ by choosing $\delta(x^{\prime})=x$ (and hence $\delta^{\prime}(x^{\prime})=(x,0)$ ) and then specializing to $x=x^{\prime}$ .

Yet, $\mathrm{ev}^{\prime}_{1}\neq\mathrm{ev}^{\prime}_{2}$ as $\mathrm{ev}^{\prime}_{1}(a)(\mathrm{swap})=(1,a)$ and $\mathrm{ev}^{\prime}_{2}(a)(\mathrm{swap})=(0,0)$ . This shows that $\mathrm{ev}^{\prime}_{1}\neq\mathrm{ev}^{\prime}_{2}$ are both “valid” semantic derivatives of the evaluation function $(\mathrm{ev})$ in our framework. In particular, it shows that semantic derivatives of higher order functions might not be unique. Our macro $\overrightarrow{\mathcal{D}}$ will return $\mathrm{ev}^{\prime}_{1}$ , but everything would still work just as well if it instead returned $\mathrm{ev}^{\prime}_{2}$ .

7.2. Canonical derivatives of higher order functions?

Differential geometers and analysts have long pursued notions of a canonical derivative of various higher order functions arising, for example, in the calculus of variations and in the study of infinite dimensional Lie groups [KM97]. Such an uncontroversial notion of derivative exists on various (infinite dimensional) spaces of functions that form suitable (so-called convenient) vector spaces, or, manifolds locally modelled on such vector spaces. At the level of generality of diffeological spaces, however, various natural notions of derivative that coincide in convenient vector spaces start to diverge and it is no longer clear what the best definition of a derivative is [CW14]. Another, fundamentally different setting that defines canonical derivatives of many higher order functions is given by synthetic differential geometry [Koc06].

While derivatives of higher order functions are of deep interest and have rightly been studied in their own right in differential geometry, we believe the situation is subtly different in computer science:

(1)

In programming applications, we use higher order programs only to construct the first order functions that we ultimately end up running and calculating derivatives of. Automatic differentiation methods can exploit this freedom: derivatives of higher order functions only matter in so far as they can be used to construct the correct derivatives of first order functions, so we can choose a simple and cheap notion of derivative among the valid options. As such, the fact that our semantics does not commit to a single notion of derivative of higher order functions can be seen as a feature rather than bug that models the pragmatics of programming practice.
(2)

While function spaces in differential geometry are typically infinite dimensional objects that are unsuitable for representation in the finite memory of a computer, higher order functions as used in programming are much more restricted: all they can do is call a function on finitely many arguments and analyse the function outputs. As such, function types in programming can be thought of as (locally) finite dimensional. In case a canonical notion of automatic derivative of higher order function is really desired, it may be worth pursuing a more intentional notion of semantics such as one based on game semantics. Such intentional techniques could capture the computational notion of higher order function better than our current (and other) extensional semantics using existing techniques from differential geometry. We hope that an exploration of such techniques might lead to an appropriate notion of computable derivative, even for higher order functions.

8. Discussion and future work

8.1. Summary

We have shown that diffeological spaces provide a denotational semantics for a higher order language with variants and inductive types (Section 4,5). We have used this to show correctness of simple forward-mode AD translations for calculating higher derivatives (Theorem 3, Theomem 8).

The structure of our elementary correctness argument for Theorem 3 is a typical logical relations proof over a denotational semantics. As explained in Section 6, this can equivalently be understood as a denotational semantics in a new kind of space obtained by categorical gluing.

Overall, then, there are two logical relations at play. One is in diffeological spaces, which ensures that all definable functions are smooth. The other is in the correctness proof (equivalently in the categorical gluing), which explicitly tracks the derivative of each function, and tracks the syntactic AD even at higher types.

8.2. Connection to the state of the art in AD implementation

As is common in denotational semantics research, we have here focused on an idealized language and simple translations to illustrate the main aspects of the method. There are a number of points where our approach is simplistic compared to the advanced current practice, as we now explain.

8.2.1. Representation of vectors

In our examples we have treated $n$ -vectors as tuples of length $n$ . This style of programming does not scale to large $n$ . A better solution would be to use array types, following [SFVPJ19]. As demonstrated by [CJS20], our categorical semantics and correctness proofs straightforwardly extend to cover them, in a similar way to our treatment of lists. In fact, [CJS20] formalizes our correctness arguments in Coq and extends them to apply to the system of [SFVPJ19].

8.2.2. Efficient forward-mode AD

For AD to be useful, it must be fast. The $(1,1)$ -AD macro $\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(1,1)}$ that we use is the basis of an efficient AD library [SFVPJ19]. Numerous optimizations are needed, ranging from algebraic manipulations, to partial evaluations, to the use of an optimizing C compiler, but the resulting implementation is performant in experiments [SFVPJ19]. The Coq formalization [CJS20] validates some of these manipulations using a similar semantics to ours. We believe the implementation in [SFVPJ19] can be extended to apply to the more general $(k,R)$ -AD methods we described in this paper through minor changes.

8.2.3. Reverse-mode and mixed-mode AD

While forward-mode AD methods are useful, many applications require reverse-mode AD, or even mixed-mode AD for efficiency. In [HSV20a], we described how our correctness proof applies to a continuation-based AD technique that closely resembles reverse-mode AD, but only has the correct complexity under a non-standard operational semantics [BMP20] (in particular, the linear factoring rule is crucial). It remains to be seen whether this technique and its correctness proof can be adapted to yield genuine reverse AD under a standard operational semantics.

Alternatively, by relying on a variation of our techniques, [Vák21] gives a correctness proof of a rather different $(1,1)$ -reverse AD algorithm that stores the (primal, adjoint)-vector pair as a struct-of-arrays rather than as an array-of-structs. Future work could explore extended its analysis to $(k,R)$ -reverse AD and mixed-mode AD for efficiently computing higher order derivatives.

8.2.4. Other language features

The idealized languages that we considered so far do not touch on several useful language constructs. For example: the use of functions that are partial (such as division) or partly-smooth (such as ReLU); phenomena such as iteration, recursion; and probabilities. Recent work by MV [Vák20] shows how our analysis of $(1,1)$ -AD extends to apply to partiality, iteration, and recursion. This development is orthogonal to the one in this paper: its methods combine directly with those in the present paper to analyze $(k,R)$ -forward mode AD of recursive programs. We leave the analysis of AD of probabilistic programs for future work.

Acknowledgment

We have benefited from discussing this work with many people, including M. Betancourt, B. Carpenter, O. Kammar, C. Mak, L. Ong, B. Pearlmutter, G. Plotkin, A. Shaikhha, J. Sigal, and others. In the course of this work, MV has also been employed at Oxford (EPSRC Project EP/M023974/1) and at Columbia in the Stan development team. This project has also received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 895827; a Royal Society University Research Fellowship; the ERC BLAST grant; the Air Force Office of Scientific Research under award number FA9550–21–1–0038; and a Facebook Research Award.

References

[AAB⁺16] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. Tensorflow: A system for large-scale machine learning. In 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pages 265–283, 2016.
[Ama12] Shun-ichi Amari. Differential-geometrical methods in statistics, volume 28. Springer Science & Business Media, 2012.
[AP20] Martín Abadi and Gordon D Plotkin. A simple differentiable programming language. In Proc. POPL 2020. ACM, 2020.
[BCLG20] Gilles Barthe, Raphaëlle Crubillé, Ugo Dal Lago, and Francesco Gavazzo. On the versatility of open logical relations: Continuity, automatic differentiation, and a containment theorem. In Proc. ESOP 2020. Springer, 2020. To appear.
[Bet18] Michael Betancourt. A geometric theory of higher-order automatic differentiation. arXiv preprint arXiv:1812.11592, 2018.
[BH11] John Baez and Alexander Hoffnung. Convenient categories of smooth spaces. Transactions of the American Mathematical Society, 363(11):5789–5825, 2011.
[BJD19] Jesse Bettencourt, Matthew J Johnson, and David Duvenaud. Taylor-mode automatic differentiation for higher-order derivatives in JAX. 2019.
[BML⁺20] Gilbert Bernstein, Michael Mara, Tzu-Mao Li, Dougal Maclaurin, and Jonathan Ragan-Kelley. Differentiating a tensor language. arXiv preprint arXiv:2008.11256, 2020.
[BMP20] Alois Brunel, Damiano Mazza, and Michele Pagani. Backpropagation in the simply typed lambda-calculus with linear negation. In Proc. POPL 2020, 2020.
[BS96] Claus Bendtsen and Ole Stauning. Fadbad, a flexible C++ package for automatic differentiation. Technical report, Technical Report IMM–REP–1996–17, Department of Mathematical Modelling, Technical University of Denmark, Lyngby, 1996.
[BS97] Claus Bendtsen and Ole Stauning. Tadiff, a flexible c++ package for automatic differentiation. TU of Denmark, Department of Mathematical Modelling, Lungby. Technical report IMM-REP-1997-07, 1997.
[CCG⁺20] J. Robin B. Cockett, Geoff S. H. Cruttwell, Jonathan Gallagher, Jean-Simon Pacaud Lemay, Benjamin MacAdam, Gordon D. Plotkin, and Dorette Pronk. Reverse derivative categories. In Proc. CSL 2020, 2020.
[CGM19] Geoff Cruttwell, Jonathan Gallagher, and Ben MacAdam. Towards formalizing and extending differential programming using tangent categories. In Proc. ACT 2019, 2019.
[CHB⁺15] Bob Carpenter, Matthew D Hoffman, Marcus Brubaker, Daniel Lee, Peter Li, and Michael Betancourt. The Stan math library: Reverse-mode automatic differentiation in C++. arXiv preprint arXiv:1509.07164, 2015.
[CJS20] Curtis Chin Jen Sem. Formalized correctness proofs of automatic differentiation in Coq. Master’s Thesis, Utrecht University, 2020. Thesis: https://dspace.library.uu.nl/handle/1874/400790. Coq code: https://github.com/crtschin/thesis.
[CRBD18] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, pages 6571–6583, 2018.
[CS96] G Constantine and T Savits. A multivariate Faa di Bruno formula with applications. Transactions of the American Mathematical Society, 348(2):503–520, 1996.
[CS11] J Robin B Cockett and Robert AG Seely. The Faa di Bruno construction. Theory and Applications of Categories, 25(15):394–425, 2011.
[CW14] J Daniel Christensen and Enxin Wu. Tangent spaces and tangent bundles for diffeological spaces. arXiv preprint arXiv:1411.5425, 2014.
[DHS11] John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12(Jul):2121–2159, 2011.
[Ell18] Conal Elliott. The simple essence of automatic differentiation. Proceedings of the ACM on Programming Languages, 2(ICFP):70, 2018.
[EM03] L Hernández Encinas and J Munoz Masque. A short proof of the generalized Faà di Bruno’s formula. Applied Mathematics Letters, 16(6):975–979, 2003.
[ER03] Thomas Ehrhard and Laurent Regnier. The differential lambda-calculus. Theoretical Computer Science, 309(1-3):1–41, 2003.
[FJL18] Roy Frostig, Matthew James Johnson, and Chris Leary. Compiling machine learning programs via high-level tracing. Systems for Machine Learning, 2018.
[FST19] Brendan Fong, David Spivak, and Rémy Tuyéras. Backprop as functor: A compositional perspective on supervised learning. In 2019 34th Annual ACM/IEEE Symposium on Logic in Computer Science (LICS), pages 1–13. IEEE, 2019.
[GSS00] Izrail Moiseevitch Gelfand, Richard A Silverman, and Richard A Silverman. Calculus of variations. Courier Corporation, 2000.
[GUW00] Andreas Griewank, Jean Utke, and Andrea Walther. Evaluating higher derivative tensors by forward propagation of univariate taylor series. Mathematics of Computation, 69(231):1117–1130, 2000.
[HG14] Matthew D Hoffman and Andrew Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1):1593–1623, 2014.
[HSV20a] Mathieu Huot, Sam Staton, and Matthijs Vákár. Correctness of automatic differentiation via diffeologies and categorical gluing. In FoSSaCS, pages 319–338, 2020.
[HSV20b] Mathieu Huot, Sam Staton, and Matthijs Vákár. Correctness of automatic differentiation via diffeologies and categorical gluing. Full version, 2020. arxiv:2001.02209.
[IZ13] Patrick Iglesias-Zemmour. Diffeology. American Mathematical Soc., 2013.
[JLS07] Peter T Johnstone, Stephen Lack, and P Sobocinski. Quasitoposes, quasiadhesive categories and Artin glueing. In Proc. CALCO 2007, 2007.
[JR11] Bart Jacobs and JMMM Rutten. An introduction to (co)algebras and (co)induction. In Advanced Topics in Bisimulation and Coinduction, pages 38–99. CUP, 2011.
[Kar01] Jerzy Karczmarczuk. Functional differentiation of computer programs. Higher-Order and Symbolic Computation, 14(1):35–57, 2001.
[KB14] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[KK04] Dana A Knoll and David E Keyes. Jacobian-free Newton–Krylov methods: a survey of approaches and applications. Journal of Computational Physics, 193(2):357–397, 2004.
[KM97] Andreas Kriegl and Peter W Michor. The convenient setting of global analysis, volume 53. American Mathematical Soc., 1997.
[Koc06] Anders Kock. Synthetic differential geometry, volume 333. Cambridge University Press, 2006.
[KSM99] Ivan Kolár, Jan Slovák, and Peter W Michor. Natural operations in differential geometry. 1999.
[KTR⁺17] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M Blei. Automatic differentiation variational inference. The Journal of Machine Learning Research, 18(1):430–474, 2017.
[KW⁺52] Jack Kiefer, Jacob Wolfowitz, et al. Stochastic estimation of the maximum of a regression function. The Annals of Mathematical Statistics, 23(3):462–466, 1952.
[Lee13] John M Lee. Smooth manifolds. In Introduction to Smooth Manifolds, pages 1–31. Springer, 2013.
[LMG18] Sören Laue, Matthias Mitterreiter, and Joachim Giesen. Computing higher order derivatives of matrix and tensor expressions. Advances in Neural Information Processing Systems, 31:2750–2759, 2018.
[LMG20] Sören Laue, Matthias Mitterreiter, and Joachim Giesen. A simple and efficient tensor calculus. In AAAI, pages 4527–4534, 2020.
[LN89] Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1-3):503–528, 1989.
[LNV21] Fernando Lucatelli Nunes and Matthijs Vákár. CHAD for expressive total languages. arXiv e-prints, pages arXiv–2110, 2021.
[LYRY20] Wonyeol Lee, Hangyeol Yu, Xavier Rival, and Hongseok Yang. On correctness of automatic differentiation for non-differentiable functions. In Advances in Neural Information Processing Systems, 2020.
[Man12] Oleksandr Manzyuk. A simply typed $\lambda$ -calculus of forward automatic differentiation. In Proc. MFPS 2012, 2012.
[Mar10] James Martens. Deep learning via Hessian-free optimization. In ICML, volume 27, pages 735–742, 2010.
[Mer04] Joel Merker. Four explicit formulas for the prolongations of an infinitesimal lie symmetry and multivariate Faa di Bruno formulas. arXiv preprint math/0411650, 2004.
[MO20] Carol Mak and Luke Ong. A differential-form pullback programming language for higher-order reverse-mode automatic differentiation. arxiv:2002.08241, 2020.
[MP21] Damiano Mazza and Michele Pagani. Automatic differentiation in PCF. Proc. ACM Program. Lang., 5(POPL):1–27, 2021. doi:10.1145/3434309.
[MS92] John C Mitchell and Andre Scedrov. Notes on sconing and relators. In International Workshop on Computer Science Logic, pages 352–378. Springer, 1992.
[Nea11] Radford M Neal. MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo, chapter 5. Chapman & Hall / CRC Press, 2011.
[PGC⁺17] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.
[Pit95] Andrew M Pitts. Categorical logic. Technical report, University of Cambridge, Computer Laboratory, 1995.
[Plo18] Gordon D Plotkin. Some principles of differential programming languages. Invited talk, POPL 2018, 2018.
[PS07] Barak A Pearlmutter and Jeffrey Mark Siskind. Lazy multivariate higher-order forward-mode ad. ACM SIGPLAN Notices, 42(1):155–160, 2007.
[PS08] Barak A Pearlmutter and Jeffrey Mark Siskind. Reverse-mode AD in a functional framework: Lambda the ultimate backpropagator. ACM Transactions on Programming Languages and Systems (TOPLAS), 30(2):7, 2008.
[Qia99] Ning Qian. On the momentum term in gradient descent learning algorithms. Neural networks, 12(1):145–151, 1999.
[RM51] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical Statistics, pages 400–407, 1951.
[Sav06] Thomas H Savits. Some statistical applications of Faa di Bruno. Journal of Multivariate Analysis, 97(10):2131–2140, 2006.
[SFVPJ19] Amir Shaikhha, Andrew Fitzgibbon, Dimitrios Vytiniotis, and Simon Peyton Jones. Efficient differentiable programming in a functional array-processing language. Proceedings of the ACM on Programming Languages, 3(ICFP):97, 2019.
[SMC20] Benjamin Sherman, Jesse Michel, and Michael Carbin. $\lambda_{S}$ : Computable semantics for differentiable programming with higher-order functions and datatypes. arXiv preprint arXiv:2007.08017, 2020.
[Sou80] Jean-Marie Souriau. Groupes différentiels. In Differential geometrical methods in mathematical physics, pages 91–128. Springer, 1980.
[Sta11] Andrew Stacey. Comparative smootheology. Theory Appl. Categ., 25(4):64–117, 2011.
[Vák20] Matthijs Vákár. Denotational correctness of foward-mode automatic differentiation for iteration and recursion. arXiv preprint arXiv:2007.05282, 2020.
[Vák21] Matthijs Vákár. Reverse AD at higher types: Pure, principled and denotationally correct. In ESOP, pages 607–634, 2021.
[VMBBL18] Bart Van Merriënboer, Olivier Breuleux, Arnaud Bergeron, and Pascal Lamblin. Automatic differentiation in ML: Where we are and where we should be going. In Advances in Neural Information Processing Systems, pages 8757–8767, 2018.
[VS21] Matthijs Vákár and Tom Smeding. CHAD: Combinatory homomorphic automatic differentiation. arXiv preprint arXiv:2103.15776, 2021.
[WGP16] Mu Wang, Assefaw Gebremedhin, and Alex Pothen. Capitalizing on live variables: new algorithms for efficient hessian computation via automatic differentiation. Mathematical Programming Computation, 8(4):393–433, 2016.
[WWE⁺19] Fei Wang, Xilun Wu, Gregory Essertel, James Decker, and Tiark Rompf. Demystifying differentiable programming: Shift/reset the penultimate backpropagator. Proceedings of the ACM on Programming Languages, 3(ICFP), 2019.
[ZHCW20] Shaopeng Zhu, Shih-Han Hung, Shouvanik Chakrabarti, and Xiaodi Wu. On the principles of differentiable quantum programming languages. In Proceedings of the 41st ACM SIGPLAN International Conference on Programming Language Design and Implementation, PLDI 2020, London, UK, June 15-20, 2020, pages 272–285. ACM, 2020. doi:10.1145/3385412.3386011.

Appendix A $\mathbf{CartSp}$ and $\mathbf{Man}$ are not cartesian closed categories

Lemma 11.

There is no continuous injection $\mathbb{R}^{d+1}\to\mathbb{R}^{d}$ .

Proof A.1.

If there were, it would restrict to a continuous injection $S^{d}\to\mathbb{R}^{d}$ . The Borsuk-Ulam theorem, however, tells us that every continuous $f:S^{d}\to\mathbb{R}^{d}$ has some $x\in S^{d}$ such that $f(x)=f(-x)$ , which is a contradiction.

Let us define the terms:

{x}_{0}:\mathbf{real},\ldots,{x}_{n}:\mathbf{real}\vdash{t}_{n}=\lambda{y}.{{x% }_{0}+{x}_{1}*y+\dots+{x}_{n}*y^{n}}:\mathbf{real}\to\mathbf{real}

Assuming that $\mathbf{CartSp}$ / $\mathbf{Man}$ is cartesian closed, observe that these get interpreted as injective continuous (because smooth) functions $\mathbb{R}^{n}\to\llbracket\mathbf{real}\to\mathbf{real}\rrbracket$ in $\mathbf{CartSp}$ and $\mathbf{Man}$ .

Theorem 12.

$\mathbf{CartSp}$ is not cartesian closed.

Proof A.2.

In case $\mathbf{CartSp}$ were cartesian closed, we would have $\llbracket\mathbf{real}\to\mathbf{real}\rrbracket=\mathbf{real}^{n}$ for some $n$ . Then, we would get, in particular a continuous injection $\llbracket{t}_{n+1}\rrbracket:\mathbb{R}^{n+1}\to\mathbb{R}^{n}$ , which contradicts Lemma 11.

Theorem 13.

$\mathbf{Man}$ is not cartesian closed.

Proof A.3.

Observe that we have $\iota_{n}:\mathbb{R}^{n}\to\mathbb{R}^{n+1}$ ; $\langle a_{0},\ldots,a_{n}\rangle\mapsto\langle a_{0},\ldots,a_{n},0\rangle$ and that $\iota_{n};\llbracket{t}_{n+1}\rrbracket=\llbracket{t}_{n}\rrbracket$ . Let us write $A_{n}$ for the image of $\llbracket{t}_{n}\rrbracket$ and $A=\cup_{n\in\mathbb{N}}A_{n}$ . Then, $A_{n}$ is connected because it is the continuous image of a connected set. Similarly, $A$ is connected because it is the non-disjoint union of connected sets. This means that $A$ lies in a single connected component of $\llbracket\mathbf{real}\to\mathbf{real}\rrbracket$ , which is a manifold with some finite dimension, say $d$ .

Take some $x\in\mathbb{R}^{d+1}$ (say, $0$ ), take some open $d$ -ball $U$ around $\llbracket{t}_{d+1}\rrbracket(x)$ , and take some open $d+1$ -ball $V$ around $x$ in $\llbracket{t}_{d+1}\rrbracket^{-1}(U)$ . Then, $\llbracket{t}_{d+1}\rrbracket$ restricts to a continuous injection from $V$ to $U$ , or equivalently, $\mathbb{R}^{d+1}$ to $\mathbb{R}^{d}$ , which contradicts Lemma 11.

	$\displaystyle Q_{\mathbf{real}}$	$\displaystyle=\{f:\mathbb{R}^{k}\to\mathbb{R}~{}\|~{}f\text{ is smooth}\}$
	$\displaystyle Q_{\boldsymbol{(}{\tau}_{1}\boldsymbol{\mathop{}}\dots% \boldsymbol{\mathop{}}{\tau}_{n}\boldsymbol{)}}$	$\displaystyle=\{\textstyle f:\mathbb{R}^{k}\to\prod_{i=1}^{n}\llfloor{\tau}_{i% }\rrceil~{}\|~{}\forall\vec{r}\in\mathbb{R}^{k}.\ \forall i.\ f_{i}(\vec{r})\in Q% _{{\tau}_{i}}\}$
	$\displaystyle Q_{{{\tau}}\to{{\sigma}}}$	$\displaystyle=\{f:\mathbb{R}^{k}\to\llfloor{\tau}\rrceil\to\llfloor{\sigma}% \rrceil~{}\|~{}\forall g\in Q_{{\tau}}.\,\lambda(\vec{r}).\,f(\vec{r})(g(\vec{r% }))\in Q_{{\sigma}}\}$

Higher Order Automatic Differentiation of Higher Order Functions

Abstract.

Key words and phrases:

1. Introduction

1.0.1. Higher order functions and differentiation.

1.0.2. From smoothness to automatic derivatives at higher types.

1.0.3. Related work and context.

1.0.4. Summary of contributions.

Relation to previous work

2. Rudiments of differentiation: how to calculate with dual numbers and Taylor approximations

2.1. First order differentiation: the chain rule and dual numbers.

Smooth functions.

2.2. Higher order differentiation: the Faà di Bruno formula and Taylor approximations.

2.3. Example: a two-dimensional second order Taylor series

2.4. Example: a one-dimensional second order Taylor series

2.5. Remark

3. A Higher Order Forward-Mode AD Translation

3.1. A simple language of smooth functions.

3.2. Syntactic automatic differentiation: a functorial macro.

Lemma 1 (Functorial macro).

Proof 3.1.

4. Semantics of differentiation

4.1. Smoothness at higher types and diffeologies

Proposition 2.

Proof 4.1.

4.2. Correctness of AD

Theorem 3 (Semantic correctness of 𝒟→→𝒟\overrightarrow{\mathcal{D}}over→ start_ARG caligraphic_D end_ARG (limited)).

Proof 4.2.

5. Extending the language: variant and inductive types

5.1. Language.

5.2. Semantics.

6. Categorical analysis of (higher order) forward AD and its correctness

6.1. Syntactic categories.

Lemma 4 (e.g. [Pit95]).

Proof 6.1 (Proof notes).

6.2. Categorical gluing and logical relations.

Proposition 5.

Proof 6.2 (Proof notes).

6.3. Correctness at all first order types, via manifolds.

Proposition 6.

Proof 6.3 (Proof notes).

Lemma 7.

Proof 6.4 (Proof notes).

Theorem 8 (Semantic correctness of 𝒟→(k,R)subscript→𝒟𝑘𝑅{\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}}over→ start_ARG caligraphic_D end_ARG start_POSTSUBSCRIPT ( italic_k , italic_R ) end_POSTSUBSCRIPT (full)).

Proof 6.5 (Proof notes).

7. Discussion: What are derivatives of higher order functions?

7.1. Automatic derivatives of higher order functions may not be unique!

Lemma 9.

Proof 7.1.

Proposition 10.

Proof 7.2.

7.2. Canonical derivatives of higher order functions?

8. Discussion and future work

8.1. Summary

8.2. Connection to the state of the art in AD implementation

8.2.1. Representation of vectors

8.2.2. Efficient forward-mode AD

8.2.3. Reverse-mode and mixed-mode AD

8.2.4. Other language features

Acknowledgment

References

Appendix A 𝐂𝐚𝐫𝐭𝐒𝐩𝐂𝐚𝐫𝐭𝐒𝐩\mathbf{CartSp}bold_CartSp and 𝐌𝐚𝐧𝐌𝐚𝐧\mathbf{Man}bold_Man are not cartesian closed categories

Lemma 11.

Proof A.1.

Theorem 12.

Proof A.2.

Theorem 13.

Proof A.3.

Theorem 3 (Semantic correctness of $\overrightarrow{\mathcal{D}}$ (limited)).

Theorem 8 (Semantic correctness of ${\scalebox{0.8}{$\overrightarrow{\mathcal{D}}$}_{(k,R)}}$ (full)).

Appendix A $\mathbf{CartSp}$ and $\mathbf{Man}$ are not cartesian closed categories