# Efficient and Sound Differentiable Programming in a Functional Array-
Processing Language

AMIR SHAIKHHA
University of Edinburgh
MATHIEU HUOT
University of Oxford
SHABNAM GHASEMIRAD
ETH Zurich
ANDREW FITZGIBBON
Graphcore
Work done while at Microsoft Research. SIMON PEYTON JONES††footnotemark:
Epic Games
DIMITRIOS VYTINIOTIS††footnotemark:
DeepMind<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Automatic differentiation (AD) is a technique for computing the derivative of
a function represented by a program. This technique is considered as the de-
facto standard for computing the differentiation in many machine learning and
optimisation software tools. Despite the practicality of this technique, the
performance of the differentiated programs, especially for functional
languages and in the presence of vectors, is suboptimal.

We present an AD system for a higher-order functional array-processing
language. The core functional language underlying this system simultaneously
supports both source-to-source forward-mode AD and global optimisations such
as loop transformations. In combination, gradient computation with forward-
mode AD can be as efficient as reverse mode, and the Jacobian matrices
required for numerical algorithms such as Gauss-Newton and Levenberg-Marquardt
can be efficiently computed.

##  1 Introduction

Functional programming (FP) and automatic differentiation (AD) have been
natural partners for sixty years, and major functional languages all have
elegant automatic differentiation packages [Elliott, 2009, Baydin et al.,
2015a, Karczmarczuk, 1999]. With the increasing importance of numerical
engineering disciplines such as machine learning, speech processing, and
computer vision, there has never been a greater need for systems which
mitigate the tedious and error-prone process of manual coding of derivatives.

Popular machine learning packages such as Tensorflow [Abadi et al., 2016] and
Pytorch [Paszke et al., 2017] greatly benefit from AD for optimisation tasks.
These systems implement embedded domain-specific languages (EDSLs) with a
predefined set of efficient combinators for manipulating tensors. Furthermore,
they can use compilation (e.g., the XLA [Leary & Wang, 2017] backend for
Tensorflow and the Glow [Rotem et al., 2018] backend for PyTorch) in order to
perform further optimisations. However, these systems are quite restrictive in
what constructs are efficiently supported; additional tensor operations are
not as efficient as the predefined set of tensor operators.

The AD systems for generic-purpose programming languages are implemented using
two modes. Forward-mode AD is relatively straightforward, both as a runtime
technique using dual numbers, or as a source-to-source program transformation.
However, forward-mode AD is usually considered wildly inefficient as a way to
compute the gradient of a function, because it involves calling the forward-
mode AD function $n$ times — and $n$ may be very large (e.g. $n=10^{6}$).

This has led to a tremendous amount of work on reverse-mode AD. As a source-
to-source transformation, reverse-mode AD is characterised by the necessity to
maintain temporary variables holding partial results, to be consumed during a
“reverse pass” of gradient computation. Modern systems devote considerable
effort (e.g., checkpointing) to optimizing the recompute/store tradeoff of
these temporaries.

Our key contribution is this: we start from the “wildly inefficient” loop in
which the forward-mode function is called $n$ times, and demonstrate that it
can be made efficient simply by applying a collection of non-AD-specific
compile-time optimising transformations. In fact, our optimisations are
sufficiently aggressive to generate code that computes the gradient with
efficiency that is sometimes _better than_ standard reverse-mode techniques.
Moreover, the entire pipeline is much simpler: there is no need to grapple
with the complexities of reverse mode — it simply falls out from the result of
optimisation.

More specifically, we take a recently introduced F# subset designed for
efficient compilation of array-processing workloads, and augment it with
vector AD primitives, yielding a functional AD tool that is competitive with
the best C/C++ tools on micro and real-world computer vision and machine
learning benchmarks, and considerably faster than others. Furthermore, we show
the correctness of the differentiation and optimisation rules by introducing
the evaluation semantics of this language.

Our contributions are as follows:

  * •

We introduce $\widetilde{\textsc{F}}$, an ML-like functional language with
array processing constructs (Section 3), and $\widetilde{\textsc{M}}$, an
embedded linear-algebra-based language.

  * •

We describe forward-mode AD (Section 2) as a simple source-to-source program
transformation for $\widetilde{\textsc{F}}$ (Section 4).

  * •

We introduce a collection of optimising transformations, none of them AD-
specific, for our language (Section 5).

  * •

We demonstrate the correctness of both source-to-source AD transformations and
optimising transformations (Section 6).

  * •

For some small but realistic benchmarks, we show that our transformations
suffice to generate extremely efficient code, starting from a naïve loop that
repeatedly calls the forward derivative of the function. We compare our
results with those of other AD systems (Section 7).

We discuss related work in Section 8.

##  2 Background

Many applications in machine learning (such as training artificial neural
networks, natural language processing, and computer vision) require computing
the derivative of an objective function. There are four techniques for
computing the derivative: 1) manual derivation of analytical derivatives, 2)
numeric differentiation, 3) symbolic differentiation, and 4) automatic
differentiation.

Among these techniques, manual derivation is error prone and time consuming,
whereas, numerical differentiation (using finite difference approximations),
although easy to implement, can lead to precision errors and is really slow
for functions with vector inputs. Symbolic differentiation involves
transforming expressions into their derivative forms. Symbolically
differentiated expressions can have an exponential overhead over the original
program if one is not careful on sharing the computation. Furthermore, this
method has no general solution for computing the derivative of programs that
include iterations and control-flow constructs [Baydin et al., 2015a,
Karczmarczuk, 1999].

Automatic Differentiation (AD) systematically applies the chain rule, and
evaluates the derivatives for the primitive arithmetic operations (such as
addition, multiplication, etc.). One of the key properties of AD is the
constant-time overhead of the differentiated program with respect to the
original code; not being careful about sharing during the process of
differentiation, can lead to code explosion [Baydin et al., 2015a].

There are two main modes for implementing AD. Forward-mode computes the
derivative part (tangent part) alongside the original computation while making
a forward pass over the program. Reverse-mode makes a forward pass to compute
the original part of the program, followed by a backward pass for computing
the derivative part (adjoint part). Consider a program in ANF [Flanagan et
al., 1993]:

l $f\big{(}x_{1},...,x_{n}\big{)}=$
let $v_{1}=e_{1}$
…
let $v_{n}=e_{n}$
$v_{n}$

To compute the derivative of this function using the forward-mode AD, we
associate a tangent variable to each variable $v_{i}$, denoted by
$\overrightarrow{v_{i}}$. Each tangent variable is computed as
$\overrightarrow{v_{i}}=\overrightarrow{x_{1}}\times\frac{\partial
v_{i}}{\partial x_{1}}+...+\overrightarrow{x_{n}}\times\frac{\partial
v_{i}}{\partial x_{n}}$. In order to compute the partial derivative of $f$
with respect to $x_{i}$, we assign $\overrightarrow{x_{i}}=1$ and
$\overrightarrow{x_{j}}=0$ for $i\neq j$, and call the transformed code.

Reverse-mode AD computes the $n$ partial derivatives simultaneously, rather
than in two different iterations. To compute the derivative of this function
using this approach, we associate an adjoint variable to each variable
$v_{i}$, denoted by $\overleftarrow{v_{i}}$, which is computed as
$\overleftarrow{v_{i}}=\frac{\partial y}{\partial v_{i}}$. As a result, if we
are interested in computing the partial derivative of function $f$ with
respect to $x_{1}$, we have to compute the value of $\overleftarrow{x_{1}}$.
Similarly, if we are interested in the partial derivative of this function
with respect to $x_{2}$, we have to compute the value of
$\overleftarrow{x_{2}}$. To do so, we have to apply the chain rule in the
reverse order.

Generally speaking, forward and reverse-mode compute a column and a row,
respectively, of the full Jacobian matrix $\mathbf{J}$ at each invocation.
More precisely, for a function with an input vector of size $n$ and an output
vector of size $m$, the forward mode approach computes a column vector of size
$m$, and the reverse mode computes a row vector of size $n$ (see Figure 1).

$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\,\,\,\hskip
14.22636pt\mathbf{J}_{f}=\frac{\partial f}{\partial x}=\leavevmode\hbox
to154.41pt{\vbox to78.6pt{\pgfpicture\makeatletter\hbox{\hskip
59.08347pt\lower-45.69896pt\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{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{}{}{{}}{{}}{{}}{
{}}{{}}\hbox{\hbox{\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{\offinterlineskip{}{}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}{{{}}}{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-26.82458pt}{-26.87119pt}\pgfsys@invoke{
}\hbox{\vbox{\halign{\pgf@matrix@init@row\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding&&\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding\cr\hfil\hskip
8.28145pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.94844pt}{0.0pt}\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{${\frac{\partial
f_{1}}{\partial x_{1}}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip
8.28145pt\hfil&\hfil\hskip 7.08301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{0.0pt}\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{${\cdots}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip
7.08301pt\hfil&\hfil\hskip 9.2253pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.89229pt}{0.0pt}\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{${\frac{\partial
f_{1}}{\partial x_{n}}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
9.2253pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip
7.08301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{0.0pt}\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{${\vdots}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
7.08301pt\hfil&\hfil\hskip 7.08301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{0.0pt}\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{${\ddots}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
7.08301pt\hfil&\hfil\hskip 7.08301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{0.0pt}\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{${\vdots}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
7.08301pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip
9.87079pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-6.53778pt}{0.0pt}\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{${\frac{\partial
f_{m}}{\partial x_{1}}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
9.87079pt\hfil&\hfil\hskip 7.08301pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.75pt}{0.0pt}\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{${\cdots}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
7.08301pt\hfil&\hfil\hskip 9.87079pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-6.53778pt}{0.0pt}\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{${\frac{\partial
f_{m}}{\partial x_{n}}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
9.87079pt\hfil\cr}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}}{{{{}}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}}\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ }}{ }
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-35.91316pt}{-27.90419pt}\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{{{\nullfont{{{
{}{}{}}}{{ {}{}{}}}}}$\left[\vbox{\hrule
height=30.40419pt,depth=30.40419pt,width=0.0pt}\right.$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}}\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{30.35759pt}{-27.90419pt}\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{{{\nullfont{{{
{}{}{}}}{{ {}{}{}}}}}$\left.\vbox{\hrule
height=30.40419pt,depth=30.40419pt,width=0.0pt}\right]$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{{}}{} {{{}}{}{
{}}{}}{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{0,0,1}\definecolor[named]{.}{rgb}{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{1}\pgfsys@invoke{
}{}\pgfsys@moveto{-25.43524pt}{20.6632pt}\pgfsys@moveto{-25.43524pt}{20.6632pt}\pgfsys@lineto{-25.43524pt}{-33.47917pt}\pgfsys@lineto{-8.47235pt}{-33.47917pt}\pgfsys@lineto{-8.47235pt}{20.6632pt}\pgfsys@closepath\pgfsys@moveto{-8.47235pt}{-33.47917pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-55.75046pt}{-42.36595pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{1}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{0,0,1}\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}Forward
Mode}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{ {}}{}{{}}{} {{{}}{}{
{}}{}}{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}{{}}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{1,0,0}\definecolor[named]{.}{rgb}{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\pgfsys@color@rgb@stroke{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{1}{0}{0}\pgfsys@invoke{
}{}\pgfsys@moveto{-25.43524pt}{20.6632pt}\pgfsys@moveto{-25.43524pt}{20.6632pt}\pgfsys@lineto{-25.43524pt}{1.4389pt}\pgfsys@lineto{26.37909pt}{1.4389pt}\pgfsys@lineto{26.37909pt}{20.6632pt}\pgfsys@closepath\pgfsys@moveto{26.37909pt}{1.4389pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{30.66417pt}{10.19403pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\pgfsys@color@rgb@stroke{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{1}{0}{0}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{1,0,0}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Reverse
Mode}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{.5,0,.5}\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}{}\pgfsys@moveto{-25.6073pt}{28.45276pt}\pgfsys@moveto{-25.6073pt}{28.45276pt}\pgfsys@lineto{-25.6073pt}{-28.45276pt}\pgfsys@lineto{-11.38092pt}{-28.45276pt}\pgfsys@lineto{-11.38092pt}{28.45276pt}\pgfsys@closepath\pgfsys@moveto{-11.38092pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{.5,0,.5}\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}{}\pgfsys@moveto{-8.5359pt}{28.45276pt}\pgfsys@moveto{-8.5359pt}{28.45276pt}\pgfsys@lineto{-8.5359pt}{-28.45276pt}\pgfsys@lineto{-5.69046pt}{-28.45276pt}\pgfsys@lineto{-5.69046pt}{28.45276pt}\pgfsys@closepath\pgfsys@moveto{-5.69046pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{.5,0,.5}\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}{}\pgfsys@moveto{-2.84544pt}{28.45276pt}\pgfsys@moveto{-2.84544pt}{28.45276pt}\pgfsys@lineto{-2.84544pt}{-28.45276pt}\pgfsys@lineto{0.0pt}{-28.45276pt}\pgfsys@lineto{0.0pt}{28.45276pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{.5,0,.5}\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}{}\pgfsys@moveto{2.84544pt}{28.45276pt}\pgfsys@moveto{2.84544pt}{28.45276pt}\pgfsys@lineto{2.84544pt}{-28.45276pt}\pgfsys@lineto{5.69046pt}{-28.45276pt}\pgfsys@lineto{5.69046pt}{28.45276pt}\pgfsys@closepath\pgfsys@moveto{5.69046pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{tikz@color}{rgb}{.5,0,.5}\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}{}\pgfsys@moveto{11.38092pt}{28.45276pt}\pgfsys@moveto{11.38092pt}{28.45276pt}\pgfsys@lineto{11.38092pt}{-28.45276pt}\pgfsys@lineto{25.6073pt}{-28.45276pt}\pgfsys@lineto{25.6073pt}{28.45276pt}\pgfsys@closepath\pgfsys@moveto{25.6073pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{21.47444pt}{-37.75476pt}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\pgfsys@color@rgb@stroke{.5}{0}{.5}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{.5}{0}{.5}\pgfsys@invoke{
}\hbox{{\definecolor[named]{.}{rgb}{.5,0,.5}\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}$\textnormal{d}\widetilde{\textsc{f}}$}}
}}\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}}$ Figure 1: The
Jacobian Matrix of a Function. Forward-mode AD computes a column of this
matrix, whereas the reverse-mode AD computes a row of this matrix.
$\textnormal{d}\widetilde{\textsc{f}}$ computes the full Jacobian matrix using
a vectorized variant of the forward-mode AD.

For a class of optimisation problems, such as various computer vision problems
using the Levenberg-Marquardt algorithm [Marquardt, 1963, Levenberg, 1944,
Moré, 1978], one is required to compute the _full_ Jacobian matrix. In such
cases, neither of the two techniques perform efficiently. To compute the full
Jacobian matrix, both forward and reverse-mode techniques must iterate over
either the columns or the rows of the Jacobian matrix, respectively. Given
that both approaches have a constant overhead over the original computation,
the forward mode technique is more efficient for computing the full Jacobian
matrix when $m\gg n$, whereas the reverse mode AD is more efficient when $n\gg
m$. However, when $n$ and $m$ are in the same range, it is not clear which
approach performs better. Moreover:

  * •

By carefully examining the body of the loops needed for computing the full
Jacobian matrix, one can observe that many computations are loop-invariant and
are unnecessarily performed multiple times. Thus, there is a lost opportunity
for loop-invariant code motion for hoisting such expressions outside the loop,
thus asymptotically improving the performance (cf. the NNMF and Bundle
Adjustment experiments in Section 7).

  * •

Furthermore, while the result of automatic differentiation is known to have
only a constant factor more arithmetic operations than the original program,
the constant can be significant; this overhead can have a dramatic impact on
the run-time performance in practice. More specifically, in applications
involving the manipulation of vectors, many intermediate vectors are allocated
that can be removed. The optimisation for eliminating such intermediate
vectors is known as deforestation [Wadler, 1988, Gill et al., 1993,
Svenningsson, 2002, Coutts et al., 2007] or loop fusion in the functional
programming community. This optimisation opens the door for many other
optimisations such as normalising loops that are iterating over sparse vectors
with a single non-zero element into a single statement (cf. Example 5 in
Section 5).

$\widetilde{\textsc{M}}$ / MATLAB / NumPyDiff. $\widetilde{\textsc{M}}$ /
MATLAB / NumPy$\widetilde{\textsc{F}}$ (Section 3.1) / F#Diff.
$\widetilde{\textsc{F}}$/ F#C / C++Diff. C / C++Lowering (Section 3.2)ADiMat
[Bischof et al., 2002] / Autograd [Maclaurin et al.,
2015]Lowering$\scriptstyle\textnormal{d}\widetilde{\textsc{f}}$ (Section 4.2)
/DiffSharp [Baydin et al., 2015b]DPS [Shaikhha et al., 2017]DPS [Shaikhha et
al., 2017]Tapenade [Hascoet & Pascual, 2013] /ADIC [Narayanan et al., 2010]
Figure 2: Compilation process in different AD systems. The solid arrows
correspond to the pipeline used in $\textnormal{d}\widetilde{\textsc{f}}$.

##  3 Overview

In this section, we start with an overview of the compilation process in
$\textnormal{d}\widetilde{\textsc{f}}$, which is shown in Figure 2. This
figure demonstrates the position of $\textnormal{d}\widetilde{\textsc{f}}$
with respect to existing AD tools.

$\textnormal{d}\widetilde{\textsc{f}}$ starts from a program written in a
high-level linear algebra library, called $\widetilde{\textsc{M}}$ (Section
3.2). This program is lowered into its implementation in a higher-order
functional language with array support, called $\widetilde{\textsc{F}}$
(Section 3.1). If a part of the program requires computing differentiation
(which are specified by using high-level differentiation API exposed by
$\textnormal{d}\widetilde{\textsc{f}}$, as mentioned in Section 4.1)
$\textnormal{d}\widetilde{\textsc{f}}$ uses AD transformation rules (Section
4.2) for transforming the involved expressions into their differentiated form.

Finally, after applying several simplifications such as loop fusion, partial
evaluation, data layout transformation, etc. (Section 5) the differentiated
program is transformed into low-level C code. The generated C code uses
efficient stack-discipline memory management by using the destination-passing
style (DPS) technique [Shaikhha et al., 2017]. Alternatively, one can use
other array programming languages such as Futhark [Henriksen et al., 2017] and
SAC [Grelck & Scholz, 2006] as the target language for differentiated
$\widetilde{\textsc{F}}$ programs.

Next, we present the core functional language used in
$\textnormal{d}\widetilde{\textsc{f}}$, on top of which we define source-to-
source AD transformation and simplification rules.

###  3.1 $\widetilde{\textsc{F}}$

$\widetilde{\textsc{F}}$ (pronounced as F smooth) is a subset of F#, an ML-
like higher-order functional programming language. It is designed to be
_expressive enough_ to make it easy to write array-processing workloads, while
simultaneously being _restricted enough_ (e.g., avoiding partially applied
functions and returning functions from lambdas which are enforced by the
(T-App) and (T-Abs) typing rules, respectively) in order to compile it to code
that is as efficient as hand-written C, with very simple and efficient memory
management [Shaikhha et al., 2017].

Syntax:

e | ::= |  e $\overline{\textnormal{{e}{}}}$ | fun $\overline{\textnormal{{x}}}$ -> e | x | – Application, Abstraction, and Variable Access
---|---|---|---
| | |  r | i | true | false | – Scalar, Index, and Boolean Values
| | |  c | – Constants (see below)
| | |  let x =e ine | – (Non-Recursive) Let Binding
| | |  ife thene elsee | – Conditional
| | |  [$\textnormal{e},...,\textnormal{e}$] | – Array Construction
T | ::= |  M | – (Non-Functional) Expression Type
| | |  $\overline{\textnormal{{T}}}$ $\Rightarrow$ M | – Function Types (No Currying)
M | ::= |  Num | – Numeric Type
| | |  Array<M> | – Vector, Matrix, … Type
| | |  M $\times$ M | – Pair Type
| | |  Bool | – Boolean Type
Num | ::= |  Double | Index | – Scalar and Index Type

Typing Rules:
(T-App) $\displaystyle\displaystyle{\hbox{\hskip
31.81651pt\vbox{\hbox{\hskip-31.81651pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}:\textnormal{$\displaystyle\textnormal{$\displaystyle\overline{\textnormal{{T}}}$}\Rightarrow\textnormal{M}$}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\overline{\textnormal{{e}}}$}:\textnormal{$\displaystyle\overline{\textnormal{{T}}}$}$}}}\vbox{}}}\over\hbox{\hskip
14.59447pt\vbox{\vbox{}\hbox{\hskip-14.59447pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\
\textnormal{$\displaystyle\overline{\textnormal{{e}}}$}:\textnormal{M}$}}}}}}$
(T-Abs) $\displaystyle\displaystyle{\hbox{\hskip
25.20836pt\vbox{\hbox{\hskip-25.20836pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\cup\textnormal{$\displaystyle\overline{\textnormal{{x}}}$}:\textnormal{$\displaystyle\overline{\textnormal{{T}}}$}\vdash\textnormal{e}{}:\textnormal{M}$}}}\vbox{}}}\over\hbox{\hskip
43.40282pt\vbox{\vbox{}\hbox{\hskip-43.40282pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{fun}}{}
{$\displaystyle\overline{\textnormal{{x}}}$} {->}
{e}}:\textnormal{$\displaystyle\textnormal{$\displaystyle\overline{\textnormal{{T}}}$}\Rightarrow\textnormal{M}$}$}}}}}}$
(T-Var) $\displaystyle\displaystyle{\hbox{\hskip
15.48611pt\vbox{\hbox{\hskip-15.48611pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{x}:\textnormal{T}\in\Gamma$}}}\vbox{}}}\over\hbox{\hskip
15.20833pt\vbox{\vbox{}\hbox{\hskip-15.20833pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{x}:\textnormal{T}$}}}}}}$
(T-Let) $\displaystyle\displaystyle{\hbox{\hskip
51.99983pt\vbox{\hbox{\hskip-51.99982pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{$\displaystyle\textnormal{{e}}_{1}$}:\textnormal{$\displaystyle\textnormal{T}{}_{1}$}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\;\;\;$}\Gamma,\textnormal{x}:\textnormal{$\displaystyle\textnormal{T}{}_{1}$}\vdash\textnormal{$\displaystyle\textnormal{{e}}_{2}$}:\textnormal{$\displaystyle\textnormal{T}{}_{2}$}$}}}\vbox{}}}\over\hbox{\hskip
47.46391pt\vbox{\vbox{}\hbox{\hskip-47.4639pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{let}}{}
x = {$\displaystyle\textnormal{{e}}_{1}$}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{in}}{}
{$\displaystyle\textnormal{{e}}_{2}$}:
{$\displaystyle\textnormal{T}{}_{2}$}}$}}}}}}$ (T-If)
$\displaystyle\displaystyle{\hbox{\hskip
50.93578pt\vbox{\hbox{\hskip-50.93578pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{1}$}:\textnormal{{Bool}}{}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\textnormal{{e}}_{2}$}:\textnormal{M}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\textnormal{{e}}_{3}$}:\textnormal{M}$}}}\vbox{}}}\over\hbox{\hskip
45.61671pt\vbox{\vbox{}\hbox{\hskip-45.6167pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{if}}{}
{$\displaystyle\textnormal{{e}}_{1}$}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{then}}{}
{$\displaystyle\textnormal{{e}}_{2}$}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{else}}{}
{$\displaystyle\textnormal{{e}}_{3}$}}:\textnormal{M}$}}}}}}$
Scalar Function Constants:
+ | - | * | / | ** : Num $\Rightarrow$ Num $\Rightarrow$ Num - : Num $\Rightarrow$ Num sin | cos | tan | log | exp : Double $\Rightarrow$ Double > | < | == | <> : Num $\Rightarrow$ Num $\Rightarrow$ Bool && | || : Bool $\Rightarrow$ Bool $\Rightarrow$ Bool not : Bool $\Rightarrow$ Bool
Vector Function Constants:
build : Index $\Rightarrow$ (Index $\Rightarrow$ M) $\Rightarrow$ Array<M>
ifold : (M $\Rightarrow$ Index $\Rightarrow$ M) $\Rightarrow$ M $\Rightarrow$
Index $\Rightarrow$ M get : Array<M>$\Rightarrow$Index $\Rightarrow$ M length
: Array<M>$\Rightarrow$Index
Pair Function Constants:
pair : $\textnormal{M}{}_{1}$ $\Rightarrow$ $\textnormal{M}{}_{2}$
$\Rightarrow$ $\textnormal{M}{}_{1}$ $\times$ $\textnormal{M}{}_{2}$ fst :
$\textnormal{M}{}_{1}$ $\times$ $\textnormal{M}{}_{2}$ $\Rightarrow$
$\textnormal{M}{}_{1}$ snd : $\textnormal{M}{}_{1}$ $\times$
$\textnormal{M}{}_{2}$ $\Rightarrow$ $\textnormal{M}{}_{2}$

Syntactic Sugar:

Matrix = Array<Array<Double>> DoubleD = Double $\times$ Double VectorD =
Array<Double $\times$ Double> MatrixD = Array<Array<Double $\times$ Double>>
$\textnormal{{e}}_{0}$[$\textnormal{{e}}_{1}$] = get $\textnormal{{e}}_{0}$
$\textnormal{{e}}_{1}$ ($\textnormal{{e}}_{0}$, $\textnormal{{e}}_{1}$) = pair
$\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$ $\textnormal{{e}}_{1}$ $bop$
$\textnormal{{e}}_{2}$ = $bop$ $\textnormal{{e}}_{1}$ $\textnormal{{e}}_{2}$
Vector = Array<Double>

Figure 3: The syntax, type system, and function constants of the core
$\widetilde{\textsc{F}}$.

Figure 3 shows the abstract syntax (parentheses can be used as necessary),
type system, and several built-in functions of $\widetilde{\textsc{F}}$.
$\overline{\textnormal{x}}$ and $\overline{\textnormal{e}}$ denote one or more
variables and expressions, respectively, which are separated by spaces,
whereas, $\overline{\textnormal{T}}$ represents one or more types which are
separated by $\Rightarrow$. In addition to the usual $\lambda$-calculus
constructs (abstraction, application, and variable access),
$\widetilde{\textsc{F}}$ supports let binding and conditionals.

$\widetilde{\textsc{F}}$ supports array programming by defining the following
built-in functions: build for producing arrays; ifold for iteration for a
particular number of times (from 0 to n-1) while maintaining a state across
iterations; length to get the size of an array; and get to index an array.

One of the key features of $\widetilde{\textsc{F}}$ is its support for both
source-to-source automatic differentiation and global optimisations such as
loop-invariant code motion and loop fusion. The transformations required for
automatic differentiation are presented in Section 4.2, and the ones for
optimisation and simplification are shown in Section 5.

Next, we show how a Linear Algebra library can be defined on top of
$\widetilde{\textsc{F}}$.

###  3.2 $\widetilde{\textsc{M}}$

$\widetilde{\textsc{M}}$ is a functional Linear Algebra library, mainly
inspired by MATLAB and R, programming languages which are heavily used by data
analysts. By providing high-level vector and matrix operations,
$\widetilde{\textsc{M}}$ frees the users from low-level details and enables
them to focus on the algorithmic aspects of the problem in hand.

Table 1: Equivalent operations in Matlab, R, NumPy, and $\widetilde{\textsc{M}}$. Matlab | R | NumPy | $\widetilde{\textsc{M}}$
---|---|---|---
A * B | A %*% B | A.dot(B) | matrixMult A B
A + B | A + B | A + B | matrixAdd A B
A’ | t(A) | A.T | matrixTranspose A
ones(n, m) | matrix(1, n, m) | ones((n, m)) | matrixOnes n m
zeros(n, m) | matrix(0, n, m) | zeros((n, m)) | matrixZeros n m
eye(n) | diag(n) | eye(n) | matrixEye n

$\widetilde{\textsc{M}}$ is simply a $\widetilde{\textsc{F}}$ library, but it
can also be thought of as an embedded domain-specific language (EDSL) [Hudak,
1996]. Figure 4 demonstrates a subset of $\widetilde{\textsc{M}}$ operations
which are defined as functions in $\widetilde{\textsc{F}}$. This library is
expressive enough for constructing vectors and matrices, element-wise
operations, accessing a slice of elements, reduction-based operations
(computing the sum of vector elements), matrix transpose, and matrix
multiplication.111We have duplicated identical implementations for functions
such as vectorMap and matrixMap (and similarly for vectorMap2 and matrixMap2)
because of the restrictions imposed by $\widetilde{\textsc{F}}$ [Shaikhha et
al., 2017]. More specifically, the C code generation process does not handle
polymorphic functions. Hence, we need to specify two different monomorphic
functions with the same implementation for such functions. Supporting more
sophisticated operations such as matrix determinant and matrix decomposition
is beyond the scope of the current paper, and we leave it for the future. As
discussed before, $\widetilde{\textsc{M}}$ is inspired by MATLAB and R. As a
result, there is a mapping among the constructs of $\widetilde{\textsc{M}}$
and these matrix-based languages. Hence, it is easily possible to translate a
program written in one of these languages to $\widetilde{\textsc{M}}$. Table 1
demonstrates the mapping among a subset of the constructs of MATLAB, R, NumPy
and $\widetilde{\textsc{M}}$.

let vectorRange = fun n ->
---
build n (fun i -> i)
let vectorFill = fun n e ->
build n (fun i -> e)
let vectorHot = fun n i ->
build n (fun j -> if i = j then 1 else 0)
let vectorMap = fun v f ->
build (length v) (fun i -> f v[i])
let vectorMap2 = fun v1 v2 f ->
build (length v1) (fun i -> f v1[i] v2[i])
let vectorZip = fun v1 v2 ->
vectorMap2 v1 v2 (pair)
let vectorAdd = fun v1 v2 ->
vectorMap2 v1 v2 (+)
let vectorEMul = fun v1 v2 ->
vectorMap2 v1 v2 ($\times$)
let vectorSMul = fun v s ->
vectorMap v (fun a -> a $\times$ s)
let vectorSum = fun v ->
ifold (fun s i -> s + v[i]) 0 (length v)
let vectorDot = fun v1 v2 ->
vectorSum (vectorEMul v1 v2)
let vectorNorm = fun v ->
sqrt (vectorDot v v)
let vectorSlice = fun v s e ->
build (e $-$ s $+$ 1) (fun i -> v[i + s])
let vectorToMatrix = fun v ->
build 1 (fun i -> v)
let vectorOutProd = fun v1 v2 ->
let m1 = vectorToMatrix v1
let m2 = vectorToMatrix v2
let m2T = matrixTranspose m2
matrixMul m1 m2T

let matrixRows = fun m -> length m
---
let matrixCols = fun m -> length (m[0])
let matrixZeros = fun r c ->
build r (fun i -> vectorFill c 0)
let matrixOnes = fun r c ->
build r (fun i -> vectorFill c 1)
let matrixEye = fun n ->
build n (fun i -> vectorHot n i)
let matrixHot = fun n m r c ->
build n (fun i ->
build m (fun j ->
if (i = r && j = c) then 1 else 0
) )
let matrixMap = fun m f ->
build (length m) (fun i -> f m[i])
let matrixMap2 = fun m1 m2 f ->
build (length m1) (fun i -> f
m1[i] m2[i])
let matrixAdd = fun m1 m2 ->
matrixMap2 m1 m2 vectorAdd
let matrixTranspose = fun m ->
build (matrixCols m) (fun i ->
build (matrixRows m) (fun j ->
m[j][i]
) )
let matrixMul = fun m1 m2 ->
let m2T = matrixTranspose m2
build (matrixRows m1) (fun i ->
build (matrixCols m2) (fun j ->
vectorDot (m1[i]) (m2T[j])
) )
let matrixTrace = fun m ->
ifold (fun s i -> s+m[i][i]) 0 (length m)

Figure 4: A subset of $\widetilde{\textsc{M}}$ constructs defined in
$\widetilde{\textsc{F}}$.

Example 1. Assume that we have a matrix $M$ and two vectors $u$ and $v$ (which
are represented as column matrices and are independent of $M$). Based on
matrix calculus one can prove that
$\frac{\partial\big{(}uMv^{T}\big{)}}{\partial M}=u^{T}v$. However, computing
the differentiated version of this function using forward-mode AD tools
requires multiple iterations over the differentiated program for every element
in the matrix $M$. By using the reverse-mode AD, one can invoke the
differentiated function only once, and the adjoint parts of the input matrix
$M$ will be filled in. We will show in Section 5 that
$\textnormal{d}\widetilde{\textsc{f}}$ derives the gradient of this expression
with respect to $M$, resulting in an expression equivalent to $u^{T}v$. This
optimises away multiple iterations over the differentiated program for each
element of matrix $M$, in contrast to the existing AD tools based on the
forward-mode AD technique.

For the moment, we only show how the matrix expression $uMv^{T}$ is expressed
in $\widetilde{\textsc{M}}$:

l let f = fun u M v ->
let um = vectorToMatrix u
let vt = matrixTranspose (vectorToMatrix v)
let m = matrixMult um (matrixMult M vt)
m[0][0]

The last expression is for accessing the single scalar element of a $1\times
1$ matrix.

$\triangle$

##  4 Differentiation

In this section, we show the differentiation process in
$\textnormal{d}\widetilde{\textsc{f}}$. First, we start by the high-level API
exposed by $\textnormal{d}\widetilde{\textsc{f}}$ to the end users. Then, we
show how $\textnormal{d}\widetilde{\textsc{f}}$ uses automatic differentiation
behind the scenes for computing derivatives. Finally, we present the
optimisations offered by $\textnormal{d}\widetilde{\textsc{f}}$, and we
demonstrate how $\textnormal{d}\widetilde{\textsc{f}}$ can use these
optimisations to deduce several matrix calculus identities.

###  4.1 High-Level API

For computing the derivative of an arbitrary function,
$\textnormal{d}\widetilde{\textsc{f}}$ provides the deriv construct. This
construct can be better thought of as a macro, which is expanded during
compilation time. The expanded expression includes the expression of the
original computation, which is given as the first argument (and can be an
arbitrary scalar, vector, or matrix expression), and the derivative of this
expression with respect to the variable given as the second argument, referred
to as the independent variable. Note that one can easily compute the
derivative of an expression with respect to a list of free variables by
multiple invocation of the deriv construct.

Figure 5 shows the implementation of the deriv construct, which is a compile-
time transformation, corresponding to computing the derivative of e with
respect to x. First, deriv constructs a lambda function which has the free
variables of the given expression as its input parameters.222deriv only
handles expressions of scalar, vector, or matrix type, and cannot be used for
function types. The same restriction applies for the extracted free variables,
but there is no restriction on the type of other sub-expressions. The produced
lambda function is given as input to source-to-source automatic
differentiation (denoted by
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$),
which can handle expressions of arbitrary type as explained later in Section
4.2. The differentiated function is applied to the dual number encoding of all
the free variables (using the
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{$i$}}\llbracket}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{}}$
compile-time transformation, corresponding to the dual number encoding of a
tensor expression of rank $i$). Based on the type of the independent variable,
the result of this applied function will be the result scalar value, or an
element of the result vector or matrix.

If the free variable is different than the input variable with respect to
which we are differentiating (i.e., the independent variable), the derivative
part is a zero scalar, vector, or matrix. Otherwise, the derivative part is a
one-hot encoding scalar, vector, or matrix. If the independent variable has a
scalar type, deriv returns the applied function. However, if the independent
variable has a vector type, deriv constructs a vector with the same number of
elements as the independent variable. For computing the $r^{\textnormal{th}}$
element of the result vector, the corresponding input vector is a one-hot
encoding with a single one at the $r^{\textnormal{th}}$ position. The
situation is similar for an independent variable with a matrix type; the
corresponding one-hot encoding matrix has a single one at the
$r^{\textnormal{th}}$ row and $c^{\textnormal{th}}$ column.

Example 2. Let us assume that we would like to compute the derivative of a
program computing the cosine function with respect to its input:

l cos a

The derivative of this program at point $a$ is represented as follows:

l snd (deriv (cos a) a)

This expression is transformed into the following expression after expanding
the deriv macro:

l snd
((${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
a -> cos
a${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(a, 1))

$\triangle$

deriv e (x: Double) =
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
v ->
e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$v${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
deriv e (x: Vector) = build (length x) (fun r ->
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
v ->
e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$v${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$)
deriv e (x: Matrix) = build (matrixRows x) (fun r -> build (matrixCols x) (fun
c ->
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
v ->
e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$v${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$))
where v = fvs(e)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$
x: Double
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
= (x, 1)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$
v: Double
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
= (v, 0)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$
v: $\overline{\textnormal{{{Double}}{}}}$
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
=
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$$\textnormal{v}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
…
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{0}}\llbracket}$$\textnormal{v}_{n}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x}}$
where v = $\textnormal{v}_{1}$ … $\textnormal{v}_{n}$
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$
x: Vector
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$
= vectorZip x (vectorHot (length x) r)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$
v: Vector
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$
= vectorZip v (vectorZeros (length v))
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$
v: $\overline{\textnormal{{{Vector}}{}}}$
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$
=
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$$\textnormal{v}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$
…
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{1}}\llbracket}$$\textnormal{v}_{n}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r}}$
where v = $\textnormal{v}_{1}$ … $\textnormal{v}_{n}$
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$
x: Matrix
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$
= matrixZip x (matrixHot (matrixRows x) (matrixCols x) r c)
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$
v: Matrix
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$
= matrixZip v (matrixZeros (matrixRows v) (matrixCols v))
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$
v: $\overline{\textnormal{{{Matrix}}{}}}$
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$
=
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$$\textnormal{v}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$
…
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{A}_{\textnormal{2}}\llbracket}$$\textnormal{v}_{n}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}_{\small\textnormal{x,r,c}}$
where v = $\textnormal{v}_{1}$ … $\textnormal{v}_{n}$

Figure 5: Implementation of the deriv construct as a source-to-source
transformation pass.

Furthermore, $\textnormal{d}\widetilde{\textsc{f}}$ provides three additional
differentiation constructs, inspired by AD tools such as DiffSharp [Baydin et
al., 2015b]: 1) diff computes the derivative of a function, from a real number
to a real number, with respect to its input, 2) grad computes the gradient of
a function, from a vector of real numbers to a real number, with respect to
its input vector, and 3) jacob computes the Jacobian matrix of a vector-valued
function, a function from a vector of real numbers to a vector of real
numbers, with respect to its input vector. Figure 6 demonstrates how these
high-level differentiation constructs are defined in terms of the source-to-
source AD transformation construct $\mathcal{D}$.

Operator | Type | Definition
---|---|---
diff |  (Double$\Rightarrow$Double) $\Rightarrow$ |  fun f x -> ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$f${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ (x, 1)
|  Double$\Rightarrow$DoubleD |
grad |  (Vector$\Rightarrow$Double) |  fun f v ->
|  $\Rightarrow$Vector$\Rightarrow$VectorD |  build (length v) (fun i ->
jacob |  (Vector$\Rightarrow$Vector) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$f${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ (vectorZip v (vectorHot (length v) i))
|  $\Rightarrow$Vector$\Rightarrow$MatrixD |  )
Figure 6: High-Level Differentiation API for $\widetilde{\textsc{F}}$.

Example 2 (Continued). For the previous example, if we would like to use the
diff construct, first we have to define the following function:

l g = fun x -> cos(x)

The derivative of this function at point $a$ is represented as follows:

l snd ((diff g) a)

which is expanded to the following program:

l snd
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$g${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(a, 1))

$\triangle$

Table 2 summarizes different matrix derivatives, and how they can be computed
using our high-level API. Note that the definition of vdiff and mdiff is
similar to diff, and the definition of mgrad is similar to grad and jacob (cf.
Figure 6). Note that the deriv construct subsumes all these operators.

Table 2: Different types of matrix derivatives. Input Type Output Type | Scalar | Vector | Matrix
---|---|---|---
Scalar | diff | vdiff | mdiff
Vector | grad | jacob | –
Matrix | mgrad | – | –

One key advantage of defining different matrix derivatives in terms of
automatic differentiation is that one no longer needs to define the matrix
calculus derivative rules for all different combinations shown in Table 2.
Instead these rules can be deduced automatically from the automatic
differentiation rules defined for scalar values. Moreover, even the algebraic
identities for matrix derivative can be deduced by using the simplification
rules presented in Section 5.

Next, we present the source code transformation required for applying
automatic differentiation rules.

###  4.2 Source-to-Source Automatic Differentiation

$\textnormal{d}\widetilde{\textsc{f}}$ relies on source-to-source translation
for implementing forward-mode automatic differentiation. Each expression is
converted into an expression containing both the original computation,
together with the derivative computation, a.k.a. the dual number technique.
The scalar expressions are transformed into a pair of values, the original
computation and the derivative computation. The vector expressions are
transformed into vectors containing tuple expressions, instead of scalar
expressions. The situation is similar for higher-rank tensors such as
matrices.

The rules for automatic differentiation are demonstrated in Figure 7.
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
specifies the AD translation for expression e. A variable y is translated to
$\dot{\textnormal{y}}$, which, as opposed to the notation used in physics,
keeps both the original and derivative part (D-Abs, D-Var, and D-Let).

The differentiation of an array of elements is an array of the differentiation
of its elements (D-Vector). Constructing an array is differentiated as an
array with the same size, however, the way that each element of the array is
constructed is differentiated (D-Build). Differentiating an iteration results
in an iteration with the same number of iterations, and with the initial state
and the next state function both differentiated (D-IFold). The differentiation
of the length and indexing an array, is the same as the length and indexing
the differentiated array, respectively (D-Length and D-Get).

For scalar expressions, we only consider the differentiation rules for real
values (i.e., expressions of type Double). For expressions of type Bool and
Index the derivative part is $\mathsf{false}$ and 0, respectively (cf. Section
6).
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
is a shorthand for extracting the original computation from the translated
term
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$,
while
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
is a shorthand for accessing the derivative part. Note that in order to avoid
redundant computation and code explosion, especially for the differentiation
rules such as the product rule (D-Mult), the arguments are bound to a new
variable.

Differentiating a pair of elements results in the pair of differentiated
elements (D-Pair). Similarly, differentiating the projection of a pair, is the
projection of the differentiated pair (D-Fst, D-Snd). For other scalar-valued
functions, the differentiation rules are similar to the corresponding rules in
mathematics.

Example 2 (Continued). In the previous example, the differentiation
transformation rules translate the program as follows:

l $\dot{\textnormal{g}}$ = fun $\dot{\textnormal{x}}$ -> -snd
($\dot{\textnormal{x}}$) * sin(fst ($\dot{\textnormal{x}}$))

Based on the definition of the diff construct, we have to use the AD version
of the function (i.e., $\dot{\textnormal{g}}$) and assign 1 to the derivative
part of the input. So the value of $\cos^{\prime}$ for the input a is computed
as follows:

l snd ((diff g) a) $\leadsto$ snd
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$g${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(a, 1)) $\leadsto$ snd ($\dot{\textnormal{g}}$ (a, 1)) $\leadsto$
-snd ((a, 1)) * sin(fst ((a, 1))) $\leadsto$ -1 * sin(a) $\leadsto$ -sin(a)

$\triangle$

Similarly, we can compute the partial derivatives of a given function, by
setting the desired derivative part to one, and the rest of derivatives to
zero. This process is illustrated in the next example.

Example 3. Assume that we would like to compute the partial derivative of the
expression a * b with respect to a, which is represented as follows in
$\widetilde{\textsc{F}}$:

l snd (deriv (a * b) a)

This expression is expanded as follows:

l snd
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
a b -> a *
b${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(a, 1) (b, 0))

Note that the derivative part of the second input is set to 0. Similar to the
previous example, the result is as follows:

l snd ((fun $\dot{\textnormal{a}}$ $\dot{\textnormal{b}}$ -> (fst
($\dot{\textnormal{a}}$)*fst ($\dot{\textnormal{b}}$), fst
($\dot{\textnormal{a}}$)*snd ($\dot{\textnormal{b}}$) + snd
($\dot{\textnormal{a}}$)*fst ($\dot{\textnormal{b}}$))) (a, 1) (b, 0))

which is evaluated as follows:

l snd ((a * b, 1 * b + a * 0)) $\leadsto$ 1 * b + a * 0 $\leadsto$ b

$\triangle$

$\lambda$-calculus Rules: |  |  |
---|---|---|---
(D-App) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Abs) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun x -> e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  fun $\dot{\textnormal{{x}}}$ -> ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-Var) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$y${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | $\dot{\textnormal{{y}}}$
(D-Let) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$let x = $\textnormal{{e}}_{1}$ in $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  let $\dot{\textnormal{{x}}}$ = ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ in ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-If) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$if $\textnormal{{e}}_{1}$ then $\textnormal{{e}}_{2}$ else $\textnormal{{e}}_{3}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  if (fst ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) then ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ else ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{3}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
Vector Rules: |  |  |
(D-Vector) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$[$\textnormal{$\textnormal{{e}}_{0}$},...,\textnormal{$\textnormal{{e}}_{n}$}$]${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | [$\textnormal{${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}${$\textnormal{{e}}_{0}$}${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$},...,\textnormal{${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}${$\textnormal{{e}}_{n}$}${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$}$]
(D-Build) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$build $\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  build (fst ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) (fun i -> (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) (i, 0))
(D-IFold) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$ifold $\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$ $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ifold (fun x i ->
|  |  | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) x (i, 0)) ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ (fst ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Get) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$[$\textnormal{{e}}_{1}$]${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)[fst ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$]
(D-Length) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$length $\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (length ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$, 0)
Pair Rules: |  |  |
(D-Pair) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$($\textnormal{{e}}_{0}$, $\textnormal{{e}}_{1}$)${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$, ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Fst) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fst $\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  fst (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Snd) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$snd $\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  snd (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$$\textnormal{{e}}_{0}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
Scalar Rules: |  |  |
|  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$, ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Add) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$ \+ $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ \+ ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-Mult) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$ * $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ \+ ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-Div) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$ / $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ \- ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) / (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$**2)
(D-Neg) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$-$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | -${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-Pow) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$ ** $\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ / ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ +
|  |  |  log (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) * ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) * (${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ ** ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Sin) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$sin $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * (cos ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Cos) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$cos $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | -${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * (sin ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
(D-Tan) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$tan $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ / ((cos ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$) ** 2)
(D-Log) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$log $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ / ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(D-Exp) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$exp $\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{E}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ * (exp ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{P}\llbracket}$$\textnormal{{e}}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$)
Typing Rules: |  |  |
(DT-Fun) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{T}{}_{1}$ $\Rightarrow$ $\textnormal{T}{}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{T}{}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ $\Rightarrow$ ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{T}{}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(DT-Num) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$Num${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = | Num $\times$ Num
(DT-Bool) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$Bool${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  Bool $\times$ Bool
(DT-Arr) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$Array<M>${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  Array<${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$M${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$>
(DT-Pair) |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{M}{}_{1}$ $\times$ $\textnormal{M}{}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ | = |  ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{M}{}_{1}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$ $\times$ ${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D_{T}}\llbracket}$$\textnormal{M}{}_{2}$${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
Figure 7: Automatic Differentiation Rules for $\widetilde{\textsc{F}}$
Expressions.

It is important to note that $\textnormal{d}\widetilde{\textsc{f}}$ performs
many of the evaluation steps shown for the previous examples during
compilation time, i.e., performs partial evaluation. Section 5 gives more
details on the optimisations and simplifications offered by
$\textnormal{d}\widetilde{\textsc{f}}$.

###  4.3 Perturbation Confusion and Nested Differentiation

In several problems such as computing the Hessian matrix, one requires to
compute the differentiation of a differentiated program. In such cases, one
should be careful dealing with tangent parts. We demonstrate this problem in
the next example.

Example 4. Here is the classical example showing the perturbation confusion
problem:

$\frac{\partial}{\partial x}\big{(}x\frac{\partial x+y}{\partial y}\big{)}$

This expression should be evaluated to $1$ at every point. However, an AD tool
can mistakenly evaluate this expression to $2$. This is because of confusing
the tangent part (perturbation) of the free variable x, while computing the
inner derivative. This is known as the perturbation confusion problem in the
AD literature. $\triangle$

If one uses the differentiation API of Figure 6, the perturbation confusion
problem appears. In order to avoid this problem, the deriv macro needs to be
used. The macro expansion of the deriv operator can be thought of as a
preprocessing step that binds each of the perturbations to a different
variable [Siskind & Pearlmutter, 2005]. Then, it is the responsibility of the
$\widetilde{\textsc{F}}$ programming language implementation (e.g., using
alpha renaming as mentioned in [Siskind & Pearlmutter, 2008]) to avoid the
perturbation confusion problem. We demonstrate this fact on our running
example.

Example 4 (Continued). The previous expression is implemented as follows in
the $\widetilde{\textsc{F}}$ language:

l fun x y ->
snd (
deriv (x * (snd (
deriv (x + y) y
))) x
)

After expanding the inner deriv macro, the following expression is derived:

l fun x y ->
snd (
deriv (
let $t_{1}$ = snd (
(fun $\dot{\textnormal{x}}$ $\dot{\textnormal{y}}$ -> (fst
($\dot{\textnormal{x}}$) + fst ($\dot{\textnormal{y}}$), snd
($\dot{\textnormal{x}}$) + snd ($\dot{\textnormal{y}}$))) (x, 0) (y, 1)
)
x * $t_{1}$) x
)

Note that the variable $t_{1}$ is created in order to avoid code explosion
that can result from the product rule (cf. Section 4.2). Expanding the outer
deriv macro results in the following expression:

l fun x y ->
snd (
(fun $\dot{\textnormal{x}}$ $\dot{\textnormal{y}}$ ->
let $\dot{\textnormal{$t_{1}$}}$ = snd (
(fun $\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$
$\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$ ->
( ( fst (fst ($\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$)) + fst (fst
($\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$)) , snd (fst
($\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$)) + snd (fst
($\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$)) ) ,
( fst (snd ($\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$)) + fst (snd
($\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$)) , snd (snd
($\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$)) + snd (snd
($\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$)) ) )
) ($\dot{\textnormal{x}}$, (0, 0)) ($\dot{\textnormal{y}}$, (1, 0))
)
(fst ($\dot{\textnormal{x}}$) * fst ($\dot{\textnormal{$t_{1}$}}$), snd
($\dot{\textnormal{x}}$) * fst ($\dot{\textnormal{$t_{1}$}}$) + fst
($\dot{\textnormal{x}}$) * snd ($\dot{\textnormal{$t_{1}$}}$))
) (x, 1) (y, 0)
)

Note that this macro expansion results in the inner perturbation variables
$\dot{\textnormal{{$\dot{\textnormal{x}}$}}}$ and
$\dot{\textnormal{{$\dot{\textnormal{y}}$}}}$ which are different from the
outer variables $\dot{\textnormal{x}}$ and $\dot{\textnormal{y}}$. This
different naming results in avoiding the perturbation confusion problem.
Finally, partially evaluating the inner expression results in the following
expression:

l fun x y -> 1

$\triangle$

##  5 Efficient Differentiation

In this section, we show how $\textnormal{d}\widetilde{\textsc{f}}$ achieves
efficient differentiable programming. First, we show several transformation
rules applicable on $\widetilde{\textsc{F}}$ expressions. We show how these
transformation rules are used to derive matrix-algebraic identities, in the
level of $\widetilde{\textsc{F}}$ expressions. Then, we show how we generate C
code from $\widetilde{\textsc{F}}$ expressions for more efficient memory
management.

###  5.1 Transformation Rules

There are various algebraic identities that one can define for
$\widetilde{\textsc{F}}$. Based on these identities, vector and matrix-based
programs, as well as differentiated programs can be heavily optimised. Figure
8 shows a set of optimisations defined for $\widetilde{\textsc{F}}$. Through
examples, we show how these rewrite rules can discover vector and matrix-level
algebraic equalities.

(fun x -> $\textnormal{{e}}_{0}$) $\textnormal{{e}}_{1}$ |  $\leadsto$ |  let x = $\textnormal{{e}}_{1}$ in $\textnormal{{e}}_{0}$
---|---|---
let x = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$ |  $\leadsto$ |  $\textnormal{{e}}_{1}$[x $\mapsto$ $\textnormal{{e}}_{0}$]
let x = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$ |  $\leadsto$ |  $\textnormal{{e}}_{1}$
(x $\not\in$ fvs($\textnormal{{e}}_{1}$))
let x = |  |  let y = $\textnormal{{e}}_{0}$ in
let y = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$ |  $\leadsto$ |  let x = $\textnormal{{e}}_{1}$
in $\textnormal{{e}}_{2}$ |  |  in $\textnormal{{e}}_{2}$
let x = $\textnormal{{e}}_{0}$ in |  |  let x = $\textnormal{{e}}_{0}$ in
let y = $\textnormal{{e}}_{0}$ in |  $\leadsto$ |  let y = x in
$\textnormal{{e}}_{1}$ |  |  $\textnormal{{e}}_{1}$
let x = $\textnormal{{e}}_{0}$ in |  |  let y = $\textnormal{{e}}_{1}$ in
let y = $\textnormal{{e}}_{1}$ in |  $\leadsto$ |  let x = $\textnormal{{e}}_{0}$ in
$\textnormal{{e}}_{2}$ |  |  $\textnormal{{e}}_{2}$
(x $\not\in$ fvs($\textnormal{{e}}_{1}$))
f(let x = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$) |  $\leadsto$ |  let x = $\textnormal{{e}}_{0}$ in f($\textnormal{{e}}_{1}$)

(a) $\lambda$-Calculus Rules

e + 0 = 0 + e |  $\leadsto$ |  e
---|---|---
e * 1 = 1 * e |  $\leadsto$ |  e
e * 0 = 0 * e |  $\leadsto$ |  0
e + -e = e - e |  $\leadsto$ |  0
$\textnormal{{e}}_{0}$ * $\textnormal{{e}}_{1}$ + $\textnormal{{e}}_{0}$ * $\textnormal{{e}}_{2}$ |  $\leadsto$ |  $\textnormal{{e}}_{0}$ * ($\textnormal{{e}}_{1}$ + $\textnormal{{e}}_{2}$)

(b) Ring-Structure Rules

(build $\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$)[$\textnormal{{e}}_{2}$] |  $\leadsto$ |  $\textnormal{{e}}_{1}$ $\textnormal{{e}}_{2}$
---|---|---
(if $\textnormal{{e}}_{0}$ $>$ $\textnormal{{e}}_{2}$)
length (build $\textnormal{{e}}_{0}$ $\textnormal{{e}}_{1}$) |  $\leadsto$ |  $\textnormal{{e}}_{0}$

(c) Loop Fusion Rules

if true |  |
---|---|---
then $\textnormal{{e}}_{1}$ |  $\leadsto$ |  $\textnormal{{e}}_{1}$
else $\textnormal{{e}}_{2}$ |  |
if false |  |
then $\textnormal{{e}}_{1}$ |  $\leadsto$ |  $\textnormal{{e}}_{2}$
else $\textnormal{{e}}_{2}$ |  |
if $\textnormal{{e}}_{0}$ |  |
then $\textnormal{{e}}_{1}$ |  $\leadsto$ |  $\textnormal{{e}}_{1}$
else $\textnormal{{e}}_{1}$ |  |
if $\textnormal{{e}}_{0}$ |  |  if $\textnormal{{e}}_{0}$
then $\textnormal{{e}}_{1}$ |  $\leadsto$ |  then $\textnormal{{e}}_{1}$[$\textnormal{{e}}_{0}$ $\mapsto$ true]
else $\textnormal{{e}}_{2}$ |  |  else $\textnormal{{e}}_{2}$[$\textnormal{{e}}_{0}$ $\mapsto$ false]
f (if $\textnormal{{e}}_{0}$ |  |  if $\textnormal{{e}}_{0}$
then $\textnormal{{e}}_{1}$ |  $\leadsto$ |  then f ($\textnormal{{e}}_{1}$)
else $\textnormal{{e}}_{2}$) |  |  else f ($\textnormal{{e}}_{2}$)

(d) Conditional Rules

fst ($\textnormal{{e}}_{0}$, $\textnormal{{e}}_{1}$) |  $\leadsto$ |  $\textnormal{{e}}_{0}$
---|---|---
snd ($\textnormal{{e}}_{0}$, $\textnormal{{e}}_{1}$) |  $\leadsto$ |  $\textnormal{{e}}_{1}$

(e) Tuple Normalisation Rules

ifold f z 0 |  $\leadsto$ |  z
---|---|---
ifold f z (n + 1) |  $\leadsto$ |  ifold (fun a i -> f a (i+1)) (f z 0) n
ifold (fun a i -> a) z n |  $\leadsto$ |  z
ifold (fun a i -> if(i = $\textnormal{{e}}_{0}$) then $\textnormal{{e}}_{1}$ else a) z n |  $\leadsto$ |  let a = z in let i = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$
(if 0 $\leq$ $\textnormal{{e}}_{0}$ $<$ n and $\textnormal{{e}}_{0}$ does not
mention a or i)

(f) Loop Normalisation Rules

ifold (fun a i -> |  |
---|---|---
($\textnormal{f}_{0}$ (fst a) i, $\textnormal{f}_{1}$ (snd a) i) |  $\leadsto$ |  (ifold $\textnormal{f}_{0}$ $\textnormal{z}_{0}$ n,ifold $\textnormal{f}_{1}$ $\textnormal{z}_{1}$ n)
) ($\textnormal{z}_{0}$, $\textnormal{z}_{1}$) n |  |

(g) Loop Fission Rule

Figure 8: Transformation Rules for $\widetilde{\textsc{F}}$. Even though none
of these rules are AD-specific, the rules of Figure 8 and Figure 8 are more
useful in the AD context.

There are various optimisations defined for scalar operations based on the
ring structure of addition and multiplication, which are shown in Figure 8.
Note that other ring-based algebraic identities, such as associativity and
commutativity, do not appear directly in the list of rules that
$\textnormal{d}\widetilde{\textsc{f}}$ applies. This is because they do not
necessarily improve the performance, unless they are combined with other
rewrite rules.

As $\widetilde{\textsc{F}}$ is based on $\lambda$-calculus, all partial
evaluation rules for this calculus come for free. Furthermore, the
optimisations defined in the literature for let-binding can also be used.
Figure 8 shows this set of rules. More specifically, the last rule subsumes
the loop-invariant code motion rule, which can have a massive impact on
performance. Consider the following C program for constructing an array:

⬇

for(int i=0; i$\langle$n; i++) {

double x = g(y);

res[i] = x * vec[i];

}

By carefully examining this program, one realises that the computation of x is
invariant to the loop, and can be hoisted outside:

⬇

double x = g(y);

for(int i=0; i$\langle$n; i++) {

res[i] = x * vec[i];

}

The corresponding $\widetilde{\textsc{F}}$ program before optimisation is as
follows:

l build n (fun i ->
let x = g y in
x * vec[i]
)

Applying the last rule of Figure 8 results in the following program:

l let x = g y in
build n (fun i -> x * vec[i])

Another source of performance issue for vector programs is the unnecessary
intermediate vectors. Consider the following program which represents adding
three vectors:

⬇

for(int i=0; i$\langle$n; i++) {

tmp[i] = vec1[i] + vec2[i];

}

for(int i=0; i$\langle$n; i++) {

res[i] = tmp[i] + vec3[i];

}

Applying loop fusion results in removing the unnecessary vector tmp, and
having a single loop instead of two:

⬇

for(int i=0; i$\langle$n; i++) {

res[i] = vec1[i] + vec2[i] + vec3[i];

}

As the vector constructs of $\widetilde{\textsc{F}}$ are based on pull arrays,
one can use the pull-array fusion rules for removing unnecessary intermediate
vectors and matrices. The two fusion rules for pull-arrays are shown in Figure
8. Apart from fusing a pipeline of vector/matrix operations, this optimisation
can also be used to derive several matrix-algebra identities.

Example 5. It is known that for a matrix $M$, the following equality holds
${(M^{T})}^{T}=M$. We show how we can derive the same equality in
$\textnormal{d}\widetilde{\textsc{f}}$. In other words, we show that:

l matrixTranspose (matrixTranspose M) = M

After let binding the inner expression, and inlining the definition of
matrixTranpose and the functions inside it, the following program is produced:

l let MT =
build (length M[0]) (fun i ->
build (length M) (fun j ->
M[j][i] ) ) in
build (length MT[0]) (fun i ->
build (length MT) (fun j ->
MT[j][i] ) )

Now, by applying the loop fusion rules (cf. Figure 8) and performing further
partial evaluation, the following expression is derived:

l build (length M) (fun i ->
build (length M[0]) (fun j ->
M[i][j] ) )

This is the same expression as M.

$\triangle$

Figure 8 shows the rewrite rules for conditional expressions. The first two
rules partially evaluate a conditional expression when its condition is
statically known. The next rule removes a conditional expression when both the
branches are the same expression. The fourth rewrite rule propagates the
result of evaluating the condition into both branches. Finally, the last
rewrite rule pushes a function applied to a conditional expression into both
branches. This results in duplicating that function, which can lead to
explosion in the size of expressions.

Figure 8 corresponds to normalisation rules for the ifold construct. The first
two rewrite rules are quite well-known; they unfold a loop the size of which
is statically known. The last two rewrite rules are more interesting and can
result in asymptotic performance improvements. The third rule turns a loop
that does not modify its state into a single statement corresponding to its
initial state. The last rule turns a loop that modifies its state only in one
of its iterations into a single statement. These two rules are especially
useful in the context of dealing with sparse vectors and matrices (e.g., one-
hot encoding vectors which can be the result of gradient computations, as well
as identity matrices) as we can see in the next example.

Example 6. It is known that for a vector $v$, the following equality holds:
$v\times I=v$, where $I$ is an identity matrix of the same dimension as $v$.
We show how $\textnormal{d}\widetilde{\textsc{f}}$ derives the same algebraic
identity. More specifically, we show that:

l let I = matrixEye (length v) in
build (length v) (fun i ->
ifold (fun a j ->
a + v[j] * I[j][i] ) 0 (length v)

is equivalent to v. Inlining and fusing this expression results in:

l build (length v) (fun i ->
ifold (fun a j ->
a + v[j] * (if (i=j) then 1 else 0) ) 0 (length v)

Applying conditional rules (cf. Figure 8) and ring-structure rules (cf. Figure
8) results in:

l build (length v) (fun i ->
ifold (fun a j ->
if (i=j) then a + v[j] else a) ) 0 (length v)

After applying the loop normalisation rules (cf. Figure 8), the following
expression is derived:

l build (length v) (fun i -> let a = 0 in let j = i in a + v[j])

Finally, performing partial evaluation and simplifications results in the
following expression:

l build (length v) (fun i -> v[i])

This is the same expression as v.

$\triangle$

Let us focus on transformations which are more related to differentiation (but
not specific to them). Many intermediate tuples resulting from the dual number
technique of AD can be removed by using partial evaluation. Figure 8 shows the
partial evaluation rules for removing the intermediate tuples which are
followed by a projection.

Partially evaluating the tuples across the boundary of a loop requires a
sophisticated analysis of the body of the loop. To simplify this task, we
perform loop fission for the loops that return a tuple of values. This is
possible only when different elements of the tuple are computed independently
in different iterations of the loop.

Consider the following C program which computes two aggregates:

⬇

for(int i=0; i$\langle$n; i++) {

res1 = res1 + vec[i];

res2 = res2 * vec[i];

}

As these two aggregates are independent of each other, loop fission can
separate their loops as follows:

⬇

for(int i=0; i$\langle$n; i++) {

res1 = res1 + vec[i];

}

for(int i=0; i$\langle$n; i++) {

res2 = res2 * vec[i];

}

Figure 8 shows how loop fission turns an iteration creating a pair of elements
into a pair of two iterations constructing independently the elements of that
pair. After performing this optimisation, if we are interested only in a
particular element of the result tuple, other loops corresponding to
irrelevant elements are removed by partial evaluation.

Based on these rewrite rules, $\textnormal{d}\widetilde{\textsc{f}}$ derives
well-known matrix calculus rules, without requiring to add a rewrite rule in
the level of matrices (i.e., $\widetilde{\textsc{M}}$). In the next section,
we suggest a strategy for the ordering of the presented rewrite rules which
supports all the examples and benchmarks reported in this article.

###  5.2 Transformation Strategy

Figure 9 shows the transformation pipeline of
$\textnormal{d}\widetilde{\textsc{f}}$. The phases shown in the first two rows
focus on the differentiated subexpression (specified by the highlighted gray
area), starting by expanding the deriv macro, followed by applying source-to-
source AD rules on the expanded expression, and then applying optimisations on
the differentiated subexpression. In the first optimisation, the body of
functions are inlined in order to open opportunities for loop fusion (Figure
8) and further partial evaluation (Figure 8). Then, we apply loop fission
(Figure 8) and tuple partial evaluation (Figure 8). Afterwards, we perform
various simplification rules consisting of algebraic identities available in
the scalar arithmetic level (e.g., associativity, commutativity, and
distributivity as shown in Figure 8), as well as $\lambda$-calculus and its
extended constructs rewrite rules (as shown in Figures 8, 8, and 8). As last
optimisation for the differentiated subexpression, we perform induction-based
loop normalisation rules of ifold (Figure 8) in order to turn loops that do
not change the state in to a single statement.

Afterwards, we focus on the whole program and again apply loop fusion and
partial evaluation in order to fuse the differentiated subexpression with the
rest of the expressions. Then, we again apply the same sequence of
optimisations as before, followed by turning the program into ANF and applying
common-subexpression elimination (CSE). The ANF makes it easier to perform
loop-invariant code motion (the last rule of Figure 8) and dead-code
elimination (DCE, which is the third rule of the same figure). Finally, after
performing simplification rules, we rely on the restrictions imposed by the
$\widetilde{\textsc{F}}$ language, such as only allowing a limited form of
recursion (using ifold and build, thus benefiting from their associated
optimisations), and produce efficient low-level C code using DPS [Shaikhha et
al., 2017].

This transformation pipeline works well in practice as we will see in Section
7. However, we do not have any guarantee that for all programs these sets of
rewrite rules with the given sequence are sufficient to make forward-mode AD
as efficient as reverse-mode AD. Also, we leave the use of alternative
transformation pipelines and using search strategies for automated rewriting
(e.g., using Monte-Carlo tree search [De Mesmay et al., 2009]) as future work.

The next example shows how $\textnormal{d}\widetilde{\textsc{f}}$ can derive a
well-known algebraic identity for the derivative of matrices using the
presented pipeline.

$\widetilde{\textsc{M}}$ Expr.MacroExpansion$\widetilde{\textsc{M}}$
Expr.${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$e${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$Differentiation$\widetilde{\textsc{M}}$
Expr.Diff. Subexpr.Loop Fusion &Partial Evaluation$\widetilde{\textsc{M}}$
Expr.Diff. Subexpr.Loop Fission &Tuple PE$\widetilde{\textsc{M}}$ Expr.Diff.
Subexpr.Simplification$\widetilde{\textsc{M}}$ Expr.Diff. Subexpr.Loop
Normalisation$\widetilde{\textsc{M}}$ Expr.Diff. Subexpr.Loop Fusion & Partial
EvaluationDiff.$\widetilde{\textsc{F}}$ Expr.Loop Fission & Tuple
PEDiff.$\widetilde{\textsc{F}}$
Expr.SimplificationDiff.$\widetilde{\textsc{F}}$
Expr.LoopNormalisationDiff.$\widetilde{\textsc{F}}$ Expr.ANF &
CSEDiff.$\widetilde{\textsc{F}}$ Expr.Code Motion&
DCEDiff.$\widetilde{\textsc{F}}$
Expr.SimplificationDiff.$\widetilde{\textsc{F}}$ Expr.DPSDiff.C Expr. Figure
9: Transformation pipeline of $\textnormal{d}\widetilde{\textsc{f}}$.

Example 7. Based on matrix calculus derivative rules, it is known that
$\frac{\partial\big{(}v_{1}\cdot v_{2}\big{)}}{\partial v_{1}}=v_{2}$, where
$\cdot$ is the vector dot product operator. We would like to show how
$\textnormal{d}\widetilde{\textsc{f}}$ can deduce the same algebraic identity.
In other words, we show that the following $\widetilde{\textsc{F}}$ program:

l vectorMap (deriv (vectorDot v1 v2) v1) snd)

is the equivalent to v2.

Macro Expansion. After expanding the deriv macro,
$\textnormal{d}\widetilde{\textsc{f}}$ produces the following program:

l vectorMap (
build (length v1) (fun i ->
${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$fun
v1 v2 -> vectorDot v1
v2${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(vectorZip v1 (vectorHot (length v1) i))
(vectorZip v2 (vectorZeros (length v2))))
) snd

Differentiation. We first focus on the highlighted expression. After applying
AD transformation rules (cf. Figure 7), the following expression is derived:

l fun $\dot{\textnormal{v1}}$ $\dot{\textnormal{v2}}$ ->
$\dot{\textnormal{vectorDot}}$ $\dot{\textnormal{v1}}$ $\dot{\textnormal{v2}}$

Loop Fusion and Partial Evaluation. Inlining the definition of
$\dot{\textnormal{vectorDot}}$ (which is derived by applying the AD
transformation rules over the $\widetilde{\textsc{F}}$ definitions given in
Figure 4) results in:

l fun $\dot{\textnormal{v1}}$ $\dot{\textnormal{v2}}$ ->
let v = build (fst (length $\dot{\textnormal{v1}}$, 0)) (fun i ->
let idx = (i, 0)
let x1 = $\dot{\textnormal{v1}}$[fst idx]
let x2 = $\dot{\textnormal{v2}}$[fst idx]
let x3 = fst x1
let x4 = snd x1
let x5 = fst x2
let x6 = snd x2
(x3*x5, x4*x5+x3*x6)
)
let z = (0, 0)
let n = (length $\dot{\textnormal{v1}}$, 0)
ifold (fun s j ->
let idx = (j, 0)
let x1 = v[fst idx]
let x2 = fst s
let x3 = snd s
let x4 = fst x1
let x5 = snd x1
(x2+x4, x3+x5)
) z (fst n)

Applying the fusion and partial evaluation rules (cf. Figure 8) results in:

l fun $\dot{\textnormal{v1}}$ $\dot{\textnormal{v2}}$ ->
ifold (fun s j ->
( (fst s) + (fst $\dot{\textnormal{v1}}$[j]) * (fst
$\dot{\textnormal{v2}}$[j]) ,
(snd s) + (snd $\dot{\textnormal{v1}}$[j]) * (fst $\dot{\textnormal{v2}}$[j])
+ (fst $\dot{\textnormal{v1}}$[j]) * (snd $\dot{\textnormal{v2}}$[j]) )
) (0, 0) (length $\dot{\textnormal{v1}}$)

Injecting the transformed expression into the initial program results in:

l vectorMap (
build (length v1) (fun i ->
(fun $\dot{\textnormal{v1}}$ $\dot{\textnormal{v2}}$ -> ifold (fun s j -> (
(fst s) + (fst $\dot{\textnormal{v1}}$[j]) * (fst $\dot{\textnormal{v2}}$[j])
, (snd s) + (snd $\dot{\textnormal{v1}}$[j]) * (fst
$\dot{\textnormal{v2}}$[j]) + (fst $\dot{\textnormal{v1}}$[j]) * (snd
$\dot{\textnormal{v2}}$[j]) ) ) (0, 0) (length $\dot{\textnormal{v1}}$))
(vectorZip v1 (vectorHot (length v1) i))
(vectorZip v2 (vectorZeros (length v2))))
) snd

Inlining the definition of vectorMap, vectorZip, vectorHot, and vectorZeros,
and applying the loop fusion, $\beta$-reduction, and partial evaluation rules
produces the following expression:

l build (length v1) (fun i ->
let s = ifold (fun s j ->
let x1 = (v1[j], if (i=j) then 1 else 0)
let x2 = (v2[j], 0)
let x3 = fst x1
let x4 = fst x2
( (fst s) + x3 * x4 ,
(snd s) + snd x1 * x4 + x3 * snd x2 )
) (0, 0) (length v1)
snd s)

After applying tuple partial evaluation (cf. Figure 8) the following program
is generated:

l build (length v1) (fun i ->
snd (ifold (fun s j ->
( (fst s) + v1[j] * v2[j] ,
(snd s) + v1[j] * 0 + (if (i=j) then 1 else 0) * v2[j] )
) (0, 0) (length v1))

Loop Fission and Tuple PE. Now we apply loop fission (cf. Figure 8):

l build (length v1) (fun i ->
snd (
ifold (fun s j -> s + v1[j] * v2[j]) 0 (length v1),
ifold (fun s j -> s + v1[j] * 0 + (if (i=j) then 1 else 0) * v2[j]) 0 (length
v1)
)
)

Note that applying the loop fission rule, does not necessarily improve the
performance; it is only after performing tuple partial evaluation rules that
the iteration responsible for the original computation is removed and the
performance is improved. Thus, the strategy for applying rewrite rules can
become tricky. For this particular rewrite rule, we only apply it, when
subsequently we can use partial evaluation to further simplify the program. To
do so, we define a compound rewrite rule that either applies these rules
together, or does not do anything. This has a similar effect to the fold-
fusion law, which can be found in the FP literature [Gibbons, 2006, Hutton,
1999]. After applying the tuple partial evaluation rule, the following program
is derived:

l build (length v1) (fun i ->
ifold (fun s j -> s + v1[j] * 0 + (if (i=j) then 1 else 0) * v2[j]) 0 (length
v1)
)

Simplification. Applying conditional rules (cf. Figure 8) and ring-based
arithmetic simplification rules (cf. Figure 8) result in:

l build (length v1) (fun i ->
(ifold (fun s j ->
if(i = j) then
(s + v2[j])
else
s) 0 (length v1)))

Loop Normalisation. By using the optimisation that turns single access
iterations into a single statement (cf. Figure 8),
$\textnormal{d}\widetilde{\textsc{f}}$ produces the following program:

l build (length v1) (fun i -> v2[i])

This program is equivalent to $v2$ if the size of the two input vectors are
the same (i.e., length $v1$ = length $v2$). Otherwise, the input program is
ill-formed.

$\triangle$

Based on the same set of transformation rules,
$\textnormal{d}\widetilde{\textsc{f}}$ derives other matrix-calculus
identities for the gradient of matrices such as $\frac{\partial
tr\big{(}M\big{)}}{\partial M}=I$, which states that the derivative of the
trace of a matrix with respect to that matrix, is an identity matrix. More
generally, $\textnormal{d}\widetilde{\textsc{f}}$ can automatically discover
the following algebraic identity if $A$ is independent of $M$: $\frac{\partial
tr\big{(}MA\big{)}}{\partial M}=A^{T}$.

Now we return to the example shown in the beginning of this paper.

Example 1 (Continued). If we have a matrix $M$ and two vectors $u$ and $v$
(which are represented as row matrices and are independent of $M$), using
matrix calculus one can prove that
$\frac{\partial\big{(}uMv^{T}\big{)}}{\partial M}=u^{T}v$. First, we start by
a partially inlined representation of this program in
$\widetilde{\textsc{F}}$:

l fun u M v ->
vectorMap (deriv (
let m =
matrixMult
(build 1 (fun i -> u))
(matrixMult M
(matrixTranspose (build 1 (fun i -> v))))
) M) snd

This expression is expanded as follows:

l let f = fun u M v ->
let m =
matrixMult
(build 1 (fun i -> u))
(matrixMult M
(matrixTranspose (build 1 (fun i -> v))))
m[0][0]
fun u M v ->
(build (length M) (fun i ->
(build (length M[0]) (fun j ->
(snd
(${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\mathcal{D}\llbracket}$f${\color[rgb]{0,0.3,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.3,0}\rrbracket}$
(vectorZip u (vectorZeros (length u)))
(matrixZip M (matrixHot (length M) (length M[0]) i j))
(vectorZip v (vectorZeros (length v)))))))))

Note that the function f is returning the only scalar element of the 1-by-1
matrix $uMv^{T}$. After performing loop fusion, loop fission and partial
evaluation the following program is derived:

l fun u M v ->
build (length M) (fun i ->
build (length M[0]) (fun j ->
u[i] * v[j]))

This program is equivalent to $u^{T}v$ if the input program is well formed,
i.e., the number of rows and columns of $M$ are the same as the length of $u$
and $v$, respectively.

$\triangle$

###  5.3 Code Generation

After applying the optimisations mentioned in the previous section, one can
further improve the efficiency by generating programs in a low-level language
with manual memory management. This way, the overhead of garbage collection
can be removed. Furthermore, by using stack-discipline memory management
techniques such as Destination-Passing Style (DPS) [Shaikhha et al., 2017],
one can benefit from efficient bump memory allocation instead of using the
expensive malloc and free calls.

Example 1 (Continued). The generated C code for the optimised differentiated
program is as follows:

⬇

matrix uMv_d(storage s, vector u, matrix M, vector v) {

matrix res = (matrix)s;

for(int r = 0; r $\langle$ M$\rightarrow\;$rows; r++) {

for(int c = 0; c $\langle$ M$\rightarrow\;$cols; c++) {

res$\rightarrow\;$elems[r][c] = u$\rightarrow\;$elems[r] *
v$\rightarrow\;$elems[c];

}

}

return res;

}

The parameter s is the storage area allocated for storing the result matrix.

$\triangle$

Up to now, we have only seen the cases where only the derivative part of the
program was of interest. If we are interested in the original part of the
program as well (e.g., the intermediate vectors cannot be fused), we need to
store both the original and derivative parts. In such cases, the
differentiated vectors, which are represented as arrays of tuples, can be
transformed into a more efficient data layout. The well-known array of structs
(AoS) to struct of arrays (SoA) transformation represents differentiated
vectors as a tuple of two numeric arrays. Further partial evaluation can
remove the unnecessary decoupled numeric arrays.

##  6 Semantics

###  6.1 Operational Semantics

We give a standard call-by-value operational semantics to the language from
Figure 3. For the purpose of this section, we turn unary functions like
$\textnormal{{sin}}:\textnormal{{Double}}{}\textnormal{$\Rightarrow$}{}\textnormal{{Double}}{}$
into typing rules, e.g.

$(\textnormal{{sin}}){\hbox{$\displaystyle\displaystyle{\hbox{\hskip
26.66669pt\vbox{\hbox{\hskip-26.66667pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{e}{}:\textnormal{{Double}}$}}}\vbox{}}}\over\hbox{\hskip
38.08334pt\vbox{\vbox{}\hbox{\hskip-38.08333pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{{sin}}(\textnormal{e}{}):\textnormal{{Double}}$}}}}}}$}}$

Similarly, for a binary function like
$\textnormal{{*}}:\textnormal{{Double}}{}\textnormal{$\Rightarrow$}{}\textnormal{{Double}}{}\textnormal{$\Rightarrow$}{}\textnormal{{Double}}{}$,
the associated typing rule is:

$(\textnormal{{*}}){\hbox{$\displaystyle\displaystyle{\hbox{\hskip
64.46649pt\vbox{\hbox{\hskip-64.46649pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{$\displaystyle\textnormal{{e}}_{1}$}:\textnormal{{Double}}\textnormal{$\displaystyle\;\;\;$}\textnormal{$\displaystyle\;\;\;$}\Gamma\vdash\textnormal{$\displaystyle\textnormal{{e}}_{2}$}:\textnormal{{Double}}$}}}\vbox{}}}\over\hbox{\hskip
37.51884pt\vbox{\vbox{}\hbox{\hskip-37.51883pt\hbox{\hbox{$\displaystyle\displaystyle\Gamma\vdash\textnormal{$\displaystyle\textnormal{{e}}_{1}$}~{}\textnormal{{*}}~{}\textnormal{$\displaystyle\textnormal{{e}}_{2}$}:\textnormal{{Double}}$}}}}}}$}}$

More generally, we denote by $\textnormal{op}_{1}$ any unary operator and by
$\textnormal{op}_{2}$ any binary one. For functions with multiple arguments,
the evaluation order is from left to right. We introduce the term $\bot$ at
every type for divergence.

Values are given as follows:

$\textnormal{v}:=\textnormal{r}~{}|~{}\textnormal{i}~{}|~{}{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{true}}{}~{}|~{}{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{false}}{}~{}|~{}\textnormal{({$\textnormal{{v}}_{0}$},
{$\textnormal{{v}}_{1}$})}~{}|~{}\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{fun}}{}
{$\overline{\textnormal{{x}}}$} {->}
{e}}~{}|~{}\texttt{[$\textnormal{v}_{0},\ldots,\textnormal{v}_{n-1}$]}$

To deal with constants like tan which are not defined everywhere, let Dom be
the domain function which assigns to a constant function its standard domain
of definition. For instance,
$\textnormal{Dom}(\textnormal{{tan}}):=\\{\textnormal{v}\in\mathbb{R}~{}|~{}\textnormal{{cos}}(\textnormal{v})\neq
0\\}$. No construct is defined on $\bot$, hence if any argument reduces to
$\bot$, the construct will also reduce to $\bot$. Constants of ground type
will be underlined in the operational semantics for clarity. For example, the
rule

$\displaystyle\displaystyle{\hbox{\hskip
27.02509pt\vbox{\hbox{\hskip-27.02507pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{v}\in\textnormal{Dom}(\textnormal{op}_{1})$}}}\vbox{}}}\over\hbox{\hskip
19.45558pt\vbox{\vbox{}\hbox{\hskip-19.45557pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{1}(\underline{\textnormal{v}})\to\underline{\textnormal{op}_{1}(\textnormal{v})}$}}}}}}$

when $\textnormal{op}_{1}=\textnormal{{tan}}$ and $\textnormal{v}=7$, leads to

$\displaystyle\displaystyle{\hbox{\hskip
28.81949pt\vbox{\hbox{\hskip-28.81949pt\hbox{\hbox{$\displaystyle\displaystyle
7\in\textnormal{Dom}(\textnormal{{tan}})$}}}\vbox{}}}\over\hbox{\hskip
20.00003pt\vbox{\vbox{}\hbox{\hskip-20.00002pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{tan}}(\underline{7})\to\underline{b}$}}}}}}$

where $b$ is the real $tan(7)$.

The reductions rules for the call-by-value small-step operational semantics
are given in Figure 10.

$\begin{array}[]{lll}{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
7.71445pt\vbox{\hbox{\hskip-7.71445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{e}\to\textnormal{e}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
18.28882pt\vbox{\vbox{}\hbox{\hskip-18.28882pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{e}{}~{}\textnormal{$\displaystyle\textnormal{{e}}_{1}$}\to\textnormal{e}^{\prime}~{}\textnormal{$\displaystyle\textnormal{{e}}_{1}$}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
10.51445pt\vbox{\hbox{\hskip-10.51445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\to\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
19.12218pt\vbox{\vbox{}\hbox{\hskip-19.12218pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{v}{}~{}\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\to\textnormal{v}{}~{}\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
44.16678pt\vbox{\vbox{}\hbox{\hskip-44.16676pt\hbox{\hbox{$\displaystyle\displaystyle(\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{fun}}{}
{$\displaystyle\overline{\textnormal{{x}}}$} {->}
{e}})\textnormal{v}\to\textnormal{e}[\textnormal{v}/\textnormal{x}]$}}}}}}$}}\end{array}$
$\begin{array}[]{lll}{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
7.71445pt\vbox{\hbox{\hskip-7.71445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{e}\to\textnormal{e}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
30.2367pt\vbox{\vbox{}\hbox{\hskip-30.2367pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{1}(\textnormal{e})\to\textnormal{op}_{1}(\textnormal{e}^{\prime})$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
10.51445pt\vbox{\hbox{\hskip-10.51445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\to\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
43.05893pt\vbox{\vbox{}\hbox{\hskip-43.05893pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{2}(\textnormal{$\displaystyle\textnormal{{e}}_{0}$},\textnormal{$\displaystyle\textnormal{{e}}_{1}$})\to\textnormal{op}_{2}(\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime},\textnormal{$\displaystyle\textnormal{{e}}_{1}$})$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
10.51445pt\vbox{\hbox{\hskip-10.51445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\to\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
41.09229pt\vbox{\vbox{}\hbox{\hskip-41.09229pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{2}(\textnormal{v},\textnormal{$\displaystyle\textnormal{{e}}_{0}$})\to\textnormal{op}_{2}(\textnormal{v},\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime})$}}}}}}$}\par}\end{array}$
$\begin{array}[]{ll}{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
10.51445pt\vbox{\hbox{\hskip-10.51445pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{$\displaystyle\textnormal{{e}}_{0}$}\to\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}$}}}\vbox{}}}\over\hbox{\hskip
64.01616pt\vbox{\vbox{}\hbox{\hskip-64.01616pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{let}}{}~{}x~{}=\textnormal{$\displaystyle\textnormal{{e}}_{0}$}~{}{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{in}}{}~{}\textnormal{$\displaystyle\textnormal{{e}}_{1}$}\to{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{let}}{}~{}x~{}=\textnormal{$\displaystyle\textnormal{{e}}_{0}$}^{\prime}~{}{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{in}}{}~{}\textnormal{$\displaystyle\textnormal{{e}}_{1}$}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
44.5974pt\vbox{\vbox{}\hbox{\hskip-44.5974pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{let}}{}~{}x~{}=\textnormal{v}~{}{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{in}}{}~{}\textnormal{e}{}\to\textnormal{e}[\textnormal{v}/x]$}}}}}}$}}\end{array}$
$\begin{array}[]{ll}{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
28.6706pt\vbox{\vbox{}\hbox{\hskip-28.6706pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{fst}}{}~{}\textnormal{({$\displaystyle\textnormal{{v}}_{0}$},
{$\displaystyle\textnormal{{v}}_{1}$})}\to\textnormal{$\displaystyle\textnormal{{v}}_{0}$}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
30.75394pt\vbox{\vbox{}\hbox{\hskip-30.75392pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{snd}}{}~{}\textnormal{({$\displaystyle\textnormal{{v}}_{0}$},
{$\displaystyle\textnormal{{v}}_{1}$})}\to\textnormal{$\displaystyle\textnormal{{v}}_{1}$}$}}}}}}$}}\end{array}$
$\begin{array}[]{ll}{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
51.04727pt\vbox{\vbox{}\hbox{\hskip-51.04727pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{if}}{}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{true}}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{then}}{}
{$\displaystyle\textnormal{{e}}_{0}$}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{else}}{}
{$\displaystyle\textnormal{{e}}_{1}$}}\to\textnormal{$\displaystyle\textnormal{{e}}_{0}$}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
51.75562pt\vbox{\vbox{}\hbox{\hskip-51.7556pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{if}}{}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{false}}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{then}}{}
{$\displaystyle\textnormal{{e}}_{0}$}
{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{else}}{}
{$\displaystyle\textnormal{{e}}_{0}$}}\to\textnormal{$\displaystyle\textnormal{{e}}_{1}$}$}}}}}}$}}\end{array}$
$\begin{array}[]{llll}{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
27.02509pt\vbox{\hbox{\hskip-27.02507pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{v}\in\textnormal{Dom}(\textnormal{op}_{1})$}}}\vbox{}}}\over\hbox{\hskip
19.45558pt\vbox{\vbox{}\hbox{\hskip-19.45557pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{1}(\underline{\textnormal{v}})\to\underline{\textnormal{op}_{1}(\textnormal{v})}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
32.16394pt\vbox{\hbox{\hskip-32.16393pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{v}\not\in\textnormal{Dom}(\textnormal{op}_{1})$}}}\vbox{}}}\over\hbox{\hskip
20.84447pt\vbox{\vbox{}\hbox{\hskip-20.84447pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{1}(\underline{\textnormal{v}})\to\bot$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
39.13063pt\vbox{\hbox{\hskip-39.13063pt\hbox{\hbox{$\displaystyle\displaystyle(\textnormal{$\displaystyle\textnormal{{v}}_{0}$},\textnormal{$\displaystyle\textnormal{{v}}_{1}$})\in\textnormal{Dom}(\textnormal{op}_{2})$}}}\vbox{}}}\over\hbox{\hskip
23.34448pt\vbox{\vbox{}\hbox{\hskip-23.34447pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{2}(\underline{\textnormal{$\displaystyle\textnormal{{v}}_{0}$}},\underline{\textnormal{$\displaystyle\textnormal{{v}}_{1}$}})\to\underline{\textnormal{op}_{2}(\textnormal{$\displaystyle\textnormal{{v}}_{0}$},\textnormal{$\displaystyle\textnormal{{v}}_{1}$})}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
42.88063pt\vbox{\hbox{\hskip-42.88063pt\hbox{\hbox{$\displaystyle\displaystyle(\textnormal{$\displaystyle\textnormal{{v}}_{0}$},\textnormal{$\displaystyle\textnormal{{v}}_{1}$})\not\in\textnormal{Dom}(\textnormal{op}_{2})$}}}\vbox{}}}\over\hbox{\hskip
24.73337pt\vbox{\vbox{}\hbox{\hskip-24.73337pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{op}_{2}(\underline{\textnormal{$\displaystyle\textnormal{{v}}_{0}$}},\underline{\textnormal{$\displaystyle\textnormal{{v}}_{1}$}})\to\bot$}}}}}}$}}\end{array}$
$\begin{array}[]{ll}{\hbox{$\displaystyle\displaystyle{\hbox{\hskip
17.44559pt\vbox{\hbox{\hskip-17.44559pt\hbox{\hbox{$\displaystyle\displaystyle
0\leq\textnormal{i}<n$}}}\vbox{}}}\over\hbox{\hskip
43.32315pt\vbox{\vbox{}\hbox{\hskip-43.32314pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{get}}}{}(\texttt{[$\displaystyle\textnormal{v}_{0},\ldots,\textnormal{v}_{n-1}$]})(\textnormal{i})\to\textnormal{v}_{i}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
38.8863pt\vbox{\vbox{}\hbox{\hskip-38.8863pt\hbox{\hbox{$\displaystyle\displaystyle\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{length}}}{\texttt{[$\displaystyle\textnormal{v}_{0},\ldots,\textnormal{v}_{n-1}$]}}\to\textnormal{n}$}}}}}}$}}\end{array}$
$\begin{array}[]{l}{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
92.03824pt\vbox{\vbox{}\hbox{\hskip-92.03824pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{build}}{}(\textnormal{n})(\textnormal{{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}{fun}}{}
{$\displaystyle\overline{\textnormal{{x}}}$} {->}
{e}})\to[e[\textnormal{0}/x],\ldots,e[\textnormal{n-1}/x]]$}}}}}}$}}\end{array}$
$\begin{array}[]{ll}{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
38.35565pt\vbox{\vbox{}\hbox{\hskip-38.35564pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{ifold}}{}(\textnormal{v})(\textnormal{$\displaystyle\textnormal{{v}}_{0}$})(\textnormal{0})\to\textnormal{$\displaystyle\textnormal{{v}}_{0}$}$}}}}}}$}}&{\hbox{$\displaystyle\displaystyle{\hbox{}\over\hbox{\hskip
80.7169pt\vbox{\vbox{}\hbox{\hskip-80.71689pt\hbox{\hbox{$\displaystyle\displaystyle{\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{ifold}}{}(\textnormal{v})(\textnormal{$\displaystyle\textnormal{{v}}_{0}$})(\textnormal{i+1})\to\textnormal{v}({\color[rgb]{0,0,0.7}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0.7}\texttt{ifold}}{}(\textnormal{v})(\textnormal{$\displaystyle\textnormal{{v}}_{0}$})(\textnormal{i}))(\textnormal{i})$}}}}}}$}}\end{array}$
Figure 10: Small-step operational semantics

Given that the operational semantics is standard, we immediately get the
following:

######  Lemma 1 (Subject reduction)

If $\Gamma\vdash e:A$ and $e\to e^{\prime}$ then $\Gamma\vdash e^{\prime}:A$.

###### Lemma 2

Every closed term of ground type $\vdash\textnormal{e}{}:\textnormal{M}{}$
reduces to a value v or diverges (i.e. reduces to $\bot$).

###  6.2 Correctness of optimisation rules

We prove correct the optimisations of Figure 8. As is usual, we denote by
$\to^{*}$ the transitive reflexive closure of $\to$.

###### Theorem 3

All the rules $\leadsto$ from figure 8 are sound. That is, if $\vdash e:A$ and
$e\leadsto e^{\prime}$, then $e\to^{*}v$ iff $e^{\prime}\to^{*}v$.

######  Proof 6.4.

By case analysis on the rule. This is straightforward for most rules so we
only present the interesting cases. $z,z_{0},z_{1}$ will reduce to values
first in all the following cases, so without loss of generality we assume that
they are values.

Rule (c).1: (build $\textnormal{{e}}_{0}$
$\textnormal{{e}}_{1}$)[$\textnormal{{e}}_{2}$]  $\leadsto$
$\textnormal{{e}}_{1}$ $\textnormal{{e}}_{2}$.
If $\textnormal{$\textnormal{{e}}_{2}$}\to^{*}\textnormal{v}$ and
$\textnormal{{e}}_{1}$$\to^{*}$ fun x -> e1’ then the RHS reduces to
$\textnormal{$\textnormal{{e}}_{1}$}^{\prime}[\textnormal{v}/x]$. If
$\textnormal{$\textnormal{{e}}_{0}$}\to^{*}\textnormal{v}^{\prime}$ then the
LHS reduces to get (build v’ (fun x -> e1’)) $\textnormal{{e}}_{2}$. It
further reduces to get [$\textnormal{{e}}_{1}$’[0/x], $\ldots$,
$\textnormal{{e}}_{1}$’[v’-1/x]] $\textnormal{{e}}_{2}$. Without loss of
generality, assume that $\textnormal{{e}}_{1}$’ is a value. Then the LHS
further reduces to get
[$\textnormal{{e}}_{1}$’[0/x],$\ldots$,$\textnormal{{e}}_{1}$’[v’-1/x]] v. By
assumption v<v’ so the arguments are in the domain of get and the LHS finally
reduces to $\textnormal{{e}}_{1}$’[v/x].

Rule (f).2: ifold f z n+1 $\leadsto$ ifold (fun a i -> f a (i+1)) (f z 0) n.
By induction on $n$. If $n=0$ then the LHS reduces to f (ifold f z 0) 0 which
reduces to f z 0. The RHS immediately reduces to f z 0. If $n>0$ then the LHS
reduces to f (ifold f z 0) n. The RHS reduces in a few steps to f ( ifold (fun
a i -> f a (i+1)) (f z 0) n-1) (n-1+1). By I.H. the first argument of f in
both cases reduces to the same value and we are done.

Rule (f).4: ifold (fun a i -> if(i = $\textnormal{{e}}_{0}$) then
$\textnormal{{e}}_{1}$ else a) z n  $\leadsto$ let a = z in let i =
$\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$ (if $\textnormal{{e}}_{0}$
does not mention a or i, and $0\leq\textnormal{$\textnormal{{e}}_{0}$}<n$).
By induction on $n>\textnormal{$\textnormal{{e}}_{0}$}$. If
$n=\textnormal{$\textnormal{{e}}_{0}$}+1$ then the LHS reduces in a few steps
to (fun a -> if n-1 = $\textnormal{{e}}_{0}$) then $\textnormal{{e}}_{1}$ else
a)(LHS(n-1)). As $n-2<\textnormal{$\textnormal{{e}}_{0}$}$ a straightforward
induction shows that LHS(n-1) reduces to z. Hence the LHS reduces to
$\textnormal{{e}}_{1}$[z/a,$\textnormal{{e}}_{0}$/i], and so does the RHS. For
$n>\textnormal{$\textnormal{{e}}_{0}$}+1$ the LHS will reduce to (fun a -> if
n-1 = $\textnormal{{e}}_{0}$) then $\textnormal{{e}}_{1}$ else a)(LHS(n-1)).
As $n>\textnormal{$\textnormal{{e}}_{0}$}+1$,
$n-1\neq\textnormal{$\textnormal{{e}}_{0}$}$ and hence the else branch will be
chosen. By I.H. LHS(n-1) reduces to the same thing as the RHS, and hence the
LHS reduces to the same thing as LHS(n-1) which reduces to the same thing as
the RHS.

Rule (g): ifold (fun a i -> ($\textnormal{f}_{0}$ (fst a) i,
$\textnormal{f}_{1}$ (snd a) i)) ($\textnormal{z}_{0}$, $\textnormal{z}_{1}$)
n  $\leadsto$ (ifold $\textnormal{f}_{0}$ $\textnormal{z}_{0}$ n, ifold
$\textnormal{f}_{1}$ $\textnormal{z}_{1}$ n).
By induction on $n$. If $n=0$, both reduce to ($\textnormal{z}_{0}$,
$\textnormal{z}_{1}$). If $n>0$, then the LHS reduces in several steps to
($\textnormal{f}_{0}$ (fst LHS(n-1)) (n-1), $\textnormal{f}_{1}$ (fst
LHS(n-1)) (n-1)), where LHS(n-1) denotes the term on the LHS for n-1 instead
of n. By I.H. ifold (fun a i -> ($\textnormal{f}_{0}$ (fst a) i,
$\textnormal{f}_{1}$ (snd a) i)) ($\textnormal{z}_{0}$, $\textnormal{z}_{1}$)
n $\to^{*}$ ($\textnormal{{v}}_{0}$, $\textnormal{{v}}_{1}$) iff the RHS at
$n-1$ reduces to ($\textnormal{{v}}_{0}$, $\textnormal{{v}}_{1}$). So the RHS
reduces to ($\textnormal{f}_{0}$ $\textnormal{{v}}_{0}$ n-1,
$\textnormal{f}_{1}$ $\textnormal{{v}}_{1}$ n-1), and so does the RHS.

Rule (a).7: f(let x = $\textnormal{{e}}_{0}$ in $\textnormal{{e}}_{1}$)
$\leadsto$ let x = $\textnormal{{e}}_{0}$ in f($\textnormal{{e}}_{1}$).