\documentclass[11pt]{article}
\usepackage{ifthen}
\usepackage{alltt}
% \usepackage{indentfirst}
\usepackage{latexsym}
\usepackage{graphicx}
\usepackage{gdefs}
% \usepackage{logdef}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{setspace}
\usepackage{rotating}
\usepackage{stmaryrd}
\usepackage{longtable}
\usepackage{hhline}
\usepackage{refs}
\usepackage{graph}
\usepackage{changebar}
\usepackage{array}
\usepackage{wrapfig}
\usepackage[small,compact]{titlesec}
\usepackage{listings}
\usepackage{verbatim}
\usepackage{hyperref}

\usepackage{tweaklist}

\renewcommand{\itemhook}{\setlength{\topsep}{0pt}%
  \setlength{\itemsep}{0pt}%
  \setlength{\parsep}{0ex}%
}

% this sets the spacings in the enumerate environment
\renewcommand{\enumhook}{\setlength{\topsep}{0pt}%
  \setlength{\itemsep}{0pt}%
  \setlength{\parsep}{0ex}%
}

% this sets the spacings in the description environment
\renewcommand{\descripthook}{\setlength{\topsep}{0pt}%
  \setlength{\itemsep}{0pt}%
  \setlength{\parsep}{0ex}%
}

\setlength{\topmargin}{0.0in}
\setlength{\headheight}{0.0in}
\setlength{\headsep}{0.0in}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0.0in}
\setlength{\evensidemargin}{0.0in}
\setlength{\textheight}{9.0in}

\setcounter{topnumber}{5}
\def\topfraction{0.99}
\def\textfraction{0.05}
\def\floatpagefraction{0.9}
\setcounter{bottomnumber}{5}
\def\bottomfraction{0.99}
\setcounter{totalnumber}{10}
\def\dbltopfraction{0.99}
\def\dblfloatpagefraction{0.8}
\setcounter{dbltopnumber}{5}

\setlength{\marginparwidth}{.95in}
\let\oldmarginpar\marginpar
\renewcommand\marginpar[1]{\-\oldmarginpar[\raggedleft\footnotesize #1]%
{\raggedright\footnotesize #1}}

\let\oldthebibliography=\thebibliography
  \let\endoldthebibliography=\endthebibliography
  \renewenvironment{thebibliography}[1]{%
    \begin{oldthebibliography}{#1}%
      \setlength{\parskip}{0ex}%
      \setlength{\itemsep}{0ex}%
  }%
  {%
    \end{oldthebibliography}%
  }


\lstset{
	basicstyle=\itshape,
	xleftmargin=3em,
	literate={->}{$\rightarrow$}{2}
	{α}{$\alpha$}{1}
	{δ}{$\delta$}{1}
}

\newcommand{\eqncenter}[1]{
{\centering 

#1

}}
\newcommand{\ITE}[3]{#1 \: ? \:  #2 : #3}
\newcommand{\ITEAligned}[3]{%
  \begin{array}{cl}%
    \multicolumn{2}{l}{{#1}} \\%
      \hspace*{2.0ex} ? & #2 \\%
      \hspace*{2.0ex} : & #3%
   \end{array}}

\newcommand{\rotxc}[1]{\begin{sideways}#1\end{sideways}}
\newcommand{\invert}[1]{\rotxc{\rotxc{#1}}}

\newcounter{LineNumber}[figure]
%\newcommand{\theHLineNumber}[0]{\arabic{figure}.\arabic{LineNumber}}
%\newcounter{LineNumber}
\newcommand{\beginLineNumber}{\setcounter{LineNumber}{0}[0]\=\+\kill}
\newcommand{\beginLLineNumber}{\setcounter{LineNumber}{0}[00]\=\+\kill}
\newcommand{\beginLLLineNumber}{\setcounter{LineNumber}{0}[000]\=\+\kill}
\newcommand{\continueLLineNumber}{[00]\=\+\kill}
\newcommand{\stopLLineNumber}{\-}
\newcommand{\LNum}{\<\refstepcounter{LineNumber}[\theLineNumber]}
\newcommand{\PLNum}[1]{\<\refstepcounter{LineNumber}[{#1}\theLineNumber]}


% Struts for use in tabular and array
\newcommand\strutA{\rule{0pt}{3.5ex}}
\newcommand\strutT{\rule{0pt}{2.6ex}}
\newcommand\strutB{\rule[-1.2ex]{0pt}{0pt}}
\newcommand\strutC{\rule[-4ex]{0pt}{0pt}}

% Frequent sets
\newcommand{\bbZ}{\mathbb{Z}}
\newcommand{\bbQ}{\mathbb{Q}}
\newcommand{\bbN}{\mathbb{N}}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbM}{\mathbb{M}}
\newcommand{\bbB}{\mathbb{B}}
\newcommand{\calE}{\mathcal{E}}
\newcommand{\calG}{\mathcal{G}}
\newcommand{\calL}{\mathcal{L}}
\newcommand{\calP}{\mathcal{P}}
\newcommand{\calR}{\mathcal{R}}
\newcommand{\calS}{\mathcal{S}}
\newcommand{\calX}{\mathcal{X}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Ss}{\mathcal{S}}
\newcommand{\SSs}{\mathcal{SS}}

% Common brackets
\newcommand{\p}[1]{\left({#1}\right)}
\newcommand{\abs}[1]{\left|{#1}\right|}
\newcommand{\norm}[1]{\left\|{#1}\right\|}
\newcommand{\sq}[1]{\left[{#1}\right]}
\newcommand{\set}[1]{\left\{ {#1} \right\}}
\newcommand{\floor}[1]{\left\lfloor{#1}\right\rfloor}
\newcommand{\ceil}[1]{\left\lceil{#1}\right\rceil}
\newcommand{\vocab}[1]{\llparenthesis #1 \rrparenthesis}

%Math symbols
\newcommand{\defeq}{\overset{\underset{\text{def}}{}}{=}}
\newcommand{\Forall}[1]{\forall {#1}\!:}
\newcommand{\Exists}[1]{\exists {#1}\!:}
\DeclareMathOperator{\lnulunion}{\curlyvee}
\DeclareMathOperator{\lnulintersect}{\curlywedge}
\DeclareMathOperator{\extend}{\otimes}
\DeclareMathOperator{\combine}{\oplus}
\DeclareMathOperator{\tensor}{\odot}
\DeclareMathOperator*{\Extend}{\bigotimes}
\DeclareMathOperator*{\Combine}{\bigoplus}
\DeclareMathOperator*{\CombinePair}{{\bigoplus_p}}
\DeclareMathOperator*{\CombineTensor}{{\textstyle{{\displaystyle \bigoplus}_{\T}}}}
\DeclareMathOperator{\Tensor}{\bigodot}
\DeclareMathOperator{\concat}{||}

\newcommand{\Readout}{\ensuremath{{\cal R}}}
\newcommand{\DeTensor}{\ensuremath{{\lightning}}}
\newcommand{\transpose}[1]{#1^t}
\newcommand{\stranspose}[1]{#1^{\OpEval{t}}}
\newcommand{\TDeTensor}{\ensuremath{{\DeTensor}^{(t,\cdot)}}}
\newcommand{\srZero}{\underline{0}}
\newcommand{\srOne}{\underline{1}}
\newcommand{\couple}{\mathcal{C}}
\newcommand{\OpEval}[1]{\llparenthesis{#1}\rrparenthesis}

% Variable names
\newcommand{\var}[1]{\ensuremath{\mathit{#1}}}
\newcommand{\proc}[1]{\ensuremath{\text{\textsc{#1}}}}
\newcommand{\startState}{\textrm{startState}}
\newcommand{\accept}{\textrm{acceptState}}
\newcommand{\reject}{\textrm{rejectState}}
\newcommand{\newState}{\textrm{newState}}
\newcommand{\Ans}{\textrm{Ans}}
\newcommand{\stateTranslationTable}{\textrm{stateTranslationTable}}
\newcommand{\callStateProxy}{\textrm{callStateProxy}}
\newcommand{\callStatesAtLevel}{\textrm{callStatesAtLevel}}
\newcommand{\callStatesAfterZero}{\textrm{callStatesAfterZero}}
\newcommand{\callStatesAfterOne}{\textrm{callStatesAfterOne}}
\newcommand{\topLevelVar}{\textrm{topLevelVar}}
\newcommand{\visited}{\textrm{visited}}
\newcommand{\even}{\textrm{even}}
\newcommand{\odd}{\textrm{odd}}
\newcommand{\Bool}{\textrm{Bool}}
\newcommand{\true}{\textrm{true}}
\newcommand{\false}{\textrm{false}}
\newcommand{\each}{\textrm{each}}
\newcommand{\internalTransSet}{\textrm{internalTransSet}}
\newcommand{\callTransSet}{\textrm{callTransSet}}
\newcommand{\returnTransSet}{\textrm{returnTransSet}}
\DeclareMathOperator{\eqeq}{\textrm{==}}
\DeclareMathOperator{\unionEqual}{\textrm{$\cup$=\,}}
\newcommand{\OBDD}{\textrm{OBDD}}
\newcommand{\unsigned}{\textrm{unsigned}}
\newcommand{\Eq}{\textrm{Eq}}
\newcommand{\Id}{\textrm{Id}}
\newcommand{\expr}{\ensuremath{\mathit{exp}}}
\newcommand{\exprt}{\ensuremath{\mathit{expt}}}
\newcommand{\Procs}{\textrm{Procs}}
\newcommand{\Coeff}{\textrm{Coeff}}
\newcommand{\Reg}{\textrm{Reg}}
\newcommand{\Rhs}{\textrm{Rhs}}
\newcommand{\Project}{\textrm{Project}}

\newcommand{\Cost}{\textit{Cost}}

\newcommand{\fONEx}{{\mathtt{f}_\mathtt{1}\mathtt{(x)}}}
\newcommand{\fTWOx}{{\mathtt{f}_\mathtt{2}\mathtt{(x)}}}
\newcommand{\fTHREEx}{{\mathtt{f}_\mathtt{3}\mathtt{(x)}}}
\newcommand{\fimmx}{{\mathtt{f}_\mathtt{i-1}\mathtt{(x)}}}
\newcommand{\fix}{{\mathtt{f}_\mathtt{i}\mathtt{(x)}}}
\newcommand{\fippx}{{\mathtt{f}_\mathtt{i+1}\mathtt{(x)}}}
\newcommand{\fjx}{{\mathtt{f}_\mathtt{j}\mathtt{(x)}}}
\newcommand{\fkx}{{\mathtt{f}_\mathtt{k}\mathtt{(x)}}}
\newcommand{\fmx}{{\mathtt{f}_\mathtt{m}\mathtt{(x)}}}
\newcommand{\fONEp}{{\mathtt{f}_\mathtt{1}'}}
\newcommand{\fTWOp}{{\mathtt{f}_\mathtt{2}'}}
\newcommand{\fTHREEp}{{\mathtt{f}_\mathtt{3}'}}
\newcommand{\fip}{{\mathtt{f}_\mathtt{i}'}}
\newcommand{\fjp}{{\mathtt{f}_\mathtt{j}'}}
\newcommand{\fmp}{{\mathtt{f}_\mathtt{m}'}}
\newcommand{\fkp}{{\mathtt{f}_\mathtt{k}'}}
\newcommand{\fONEpx}{{\mathtt{f}_\mathtt{1}'\mathtt{(x)}}}
\newcommand{\fTWOpx}{{\mathtt{f}_\mathtt{2}'\mathtt{(x)}}}
\newcommand{\fTHREEpx}{{\mathtt{f}_\mathtt{3}'\mathtt{(x)}}}
\newcommand{\fipx}{{\mathtt{f}_\mathtt{i}'\mathtt{(x)}}}


\title{Automatic Differentiation and Backpropagation \\
       CS701 \\
}

\author{Thomas Reps \\ \small[Based, in part, on notes taken by Naveen Neelakandan on December 1, 2015]}
\date{}

\begin{document}

\maketitle
\begin{abstract}
This lecture discusses the relationship between \emph{automatic
differentiation} and \emph{backpropagation}.
Automatic differentiation (AD) is a technique that takes an
implementation of a numerical function $f$ (computed using
floating-point numbers) and creates an implementation of $f'$.
We explain several techniques for performing AD.
For forward-mode AD, we give an explicit transformation
of the program, as well as a way to implement AD using
operator overloading.
We then reformulate AD as a path problem over a so-called
\emph{computation graph}.
The path problem can be solved in either direction:
one direction corresponds to forward-mode AD;
the other corresponds to reverse-mode AD.

Finally, we show how reverse-mode AD forms the basis of the
backpropagation algorithm for training a neural net:
in backpropagation, the weight adjustments computed for a neural net,
a given input, and an expected output are found by performing
reverse-mode AD on the computation graph for (i) the neural net,
along with (ii) additional computation-graph elements for the
sum-of-squared differences of the neural net's output with the expected output.
\end{abstract}

\section{Introduction}
\label{Se:Introduction}

This lecture provides another payoff for our having invested so much time
on path problems.
The point of the lecture is to show you how the core step used to
train a neural net---namely, \emph{backpropagation}---involves solving a path problem.
Our path to presenting this concept is a bit indirect, but hopefully enlightening:
\begin{itemize}
  \item
    forward-mode automatic differentiation
  \item
    reverse-mode automatic differentiation
  \item
    expressing a neural net as an expression DAG
  \item
    backpropagation
\end{itemize}

\section{Forward-Mode Automatic Differentiation}
\label{Se:ForwardModeAutomaticDifferentiation}

This section\footnote{This section is adapted from \cite[\S2]{HOSC:RR03}.}
describes so-called \emph{forward-mode automatic differentiation},
which provides a way to compute, with high accuracy, the derivative
of a function defined by a program.
An AD tool transforms a program that computes a numerical function $F(x)$
into a related program that computes the derivative $F'(x)$.
These tools address the following issue:
Suppose that you have a program {\tt F(x)\/} that computes a numerical
function $F(x)$.
It is a very bad idea to try to compute $F'(x_0)$, the
value of the derivative of $F$ at $x_0$, by picking a small value
{\tt delta\_x\/} and invoking the following program with the argument
$\mathtt{x}_\mathtt{0}$:\footnote{In the remainder of this section,
  {\tt Courier Font\/} is used to denote
  functions defined by programs, whereas {\em Italic Font\/} is used to
  denote mathematical functions.
  That is, $F(x)$ denotes a function (evaluated over real numbers),
  whereas {\tt F(x)\/} denotes a program (evaluated over floating-point
  numbers).
  We adhere to this convention both in concrete examples that involve
  C++ code, as well as in more abstract discussions in order to
  distinguish between a mathematical function and a program that
  implements the function.

  The example programs are all written in C++,
  although the ideas described apply to other programming
  languages---including functional programming languages
  (cf. \cite{icfp:Karc99,HOSC:Karc01})---as well as to other
  imperative languages.
  To emphasize the links between mathematical concepts and their
  implementations in C++, we take the liberty of sometimes using
  $'$ and/or subscripts on C++ identifiers.}
\begin{equation}
  \label{Box:FPrimeNaive}
  \text{
  \begin{tabular}{|l|}
    \cline{1-1}
        {\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          float delta\_x = \ldots $\B{\mbox{some small value}}$ \ldots ; \\
          fl\=oat F$'$\_{naive}(float x)\{ \+ \\
            return (F(x + delta\_x) - F(x))/delta\_x; \- \\
          \}
        \end{tabbing}
        \end{minipage}
        } \\
    \cline{1-1}
  \end{tabular}
  }
\end{equation}
For a small enough value of {\tt delta\_x\/}, the values of
{\tt F($\mathtt{x}_\mathtt{0}$+delta\_x)\/} and
{\tt F($\mathtt{x}_\mathtt{0}$)\/} will usually be very close.
Round-off errors in the computation of
{\tt F($\mathtt{x}_\mathtt{0}$+delta\_x)\/} and
{\tt F($\mathtt{x}_\mathtt{0}$)\/} are magnified
by the subtraction of the two quantities, and further
amplified by the division by the small quantity
{\tt delta\_x\/}, which may cause the overall result
to be useless.
Automatic differentiation sidesteps this problem
by computing derivatives in another fashion.

Automatic differentiation can be illustrated by means
of the following example:

\begin{Exa}\label{Exa:ProdExample}\cite{Zippel-7-96-PC-TWR}.
Suppose that we have been given a collection of programs
$\mathtt{f}_\mathtt{i}$ for the functions $f_i$, $1 \leq i \leq k$,
together with the program {\tt Prod\/} shown below,
which computes the function $Prod(x) = \prod_{i=1}^k {f_i}(x)$.
In addition, suppose that we have also been given programs
$\fip$ for the functions $f_i'$, $1 \leq i \leq k$.
Finally, suppose that we wish to obtain a program $\mathtt{Prod}'$
that computes the function $Prod'(x)$.
Column two of the table given below shows mathematical expressions for
$Prod(x)$ and $Prod'(x)$.
Column three shows two C++ procedures: Procedure {\tt Prod\/} computes $Prod(x)$;
procedure $\mathtt{Prod}'$ is the procedure that an
automatic-differentiation system would create to compute
$Prod'(x)$.
\begin{center}
%\resizebox{\textwidth}{!}{
  \begin{tabular}{|@{\hspace{.3ex}}c@{\hspace{.3ex}}|@{\hspace{.3ex}}c@{\hspace{.3ex}}|@{\hspace{.3ex}}c@{\hspace{0ex}}|}
     \cline{1-3}
               & Mathematical Notation & Programming Notation \\
     \cline{1-3}
        Function
      & $Prod(x) = {\displaystyle \prod_{i=1}^k } {f_i}(x)$
      & {\footnotesize{\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          fl\=oat Prod(float x)\{ \+ \\
            float ans = 1.0; \\
            fo\=r (int i = 1; i <= k; i++)\{ \+ \\
              ans = ans * $\fix$; \- \\
            \} \\
            return ans; \- \\
          \}
        \end{tabbing}
        \end{minipage}
        }} \\
     \cline{1-3}
        Derivative
      & $Prod'(x) = {\displaystyle \sum_{i=1}^k } f_i'(x) * {\displaystyle \prod_{j \neq i} } f_j(x)$
      & {\footnotesize{\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          fl\=oat Prod$'$(float x)\{ \+ \\
            float ans$'$ = 0.0; \\
            float ans = 1.0; \\
            fo\=r (int i = 1; i <= k; i++)\{ \+ \\
              ans$'$ = ans$'$ * $\fix$ + ans * $\fipx$; \\
              ans = ans * $\fix$; \- \\
            \} \\
            return ans$'$; \- \\
          \}
        \end{tabbing}
        \end{minipage}
        }} \\
     \cline{1-3}
  \end{tabular}
%}
\end{center}
Notice that program $\mathtt{Prod}'$ resembles program {\tt Prod\/},
as opposed to $\mathtt{F}'\mathtt{\_naive}$ (see box~(\ref{Box:FPrimeNaive})).
$\mathtt{Prod}'$ preserves accuracy in its computation of the
derivative because, as illustrated below in
Example~\ref{Exa:ProdTraceExample}, it is based on the rules for the
exact computation of derivatives, rather than on the kind of
computation performed by $\mathtt{F}'\mathtt{\_naive}$.
\end{Exa}

The transformation illustrated above is merely one instance of a
general transformation that can be applied to any program:
Given a program {\tt G\/} as input, the transformation produces a
derivative-computing program $\mathtt{G}'$.
The method for constructing $\mathtt{G}'$ is as follows:
\begin{itemize}
  \item
    For each variable {\tt v\/} of type {\tt float\/} used in {\tt G\/},
    another {\tt float\/} variable $\mathtt{v}'$ is introduced.
  \item
    Each statement in {\tt G\/} of the form
    ``{\tt v = $exp$;\/}'',
    where $exp$ is an arithmetic expression, is transformed into
    ``{\tt v$'$ = $exp'$; v = $exp$;\/}'',
    where $exp'$ is the expression for the derivative of $exp$.
    If $exp$ involves calls to a procedure {\tt g\/}, then $exp'$
    may involve calls to both {\tt g\/} and $\mathtt{g}'$.
  \item
    Each return statement in {\tt G\/} of the form
    ``{\tt return v;\/}'' is transformed into
    ``{\tt return v$'$;\/}''.
\end{itemize}
In general, this transformation can be justified by appealing
to the chain rule of differential calculus (see below).

\begin{Exa}
\label{Exa:ProdTraceExample}
For Example~\ref{Exa:ProdExample}, we can demonstrate the
correctness of the transformation by symbolically executing
{\tt Prod\/}$'$ for a few iterations, comparing the values of
{\tt ans\/}$'$ and {\tt ans\/} (as functions of {\tt x\/}) at the
start of each iteration of the for-loop:
\begin{center}
%\resizebox{\textwidth}{!}{
  \begin{tabular}{|c|c|c|}
    \cline{1-3}
      Iteration & Value of $\mathtt{ans}'$ (as a function of {\tt x\/}) & Value of {\tt ans\/} \\
                &                                                       & (as a function of {\tt x\/}) \\
    \cline{1-3}
      0 & {\tt 0.0\/} & {\tt 1.0\/} \\
    \cline{1-3}
      1 & $\fONEpx$ & $\fONEx$ \\
    \cline{1-3}
      2 & $\fONEpx * \fTWOx + \fONEx * \fTWOpx$ & $\fONEx * \fTWOx$ \\
    \cline{1-3}
      3 & $\begin{array}{rl}
               & \fONEpx * \fTWOx * \fTHREEx \\
             + & \fONEx * \fTWOpx * \fTHREEx \\
             + & \fONEx * \fTWOx * \fTHREEpx
          \end{array}$
        & $\fONEx * \fTWOx * \fTHREEx$ \\
    \cline{1-3}
      $\ldots$ & $\ldots$ & $\ldots$ \\
    \cline{1-3}
      $k$ & ${\displaystyle \sum_{i=1}^k } \fipx * {\displaystyle \prod_{j \neq i} } \fjx$ & ${\displaystyle \prod_{i=1}^k } \fix$ \\
    \cline{1-3}
  \end{tabular}
%}
\end{center}
The loop maintains the invariant that, at the start of
each iteration,
$\mathtt{ans}'\mathtt{(x)} = \frac{d}{d\mathtt{x}} \mathtt{ans(x)}$.

The value of $\mathtt{ans}'$ on the $3^{rd}$
iteration would actually be computed with the terms grouped as follows:
\begin{equation}
  \label{Eq:ProdTraceExample}
  \mathtt{(} \fONEpx \mathtt{*} \fTWOx  \mathtt{+}  \fONEx \mathtt{*} \fTWOpx \mathtt{)} \mathtt{*} \fTHREEx \mathtt{+} \mathtt{(} \fONEx \mathtt{*} \fTWOx \mathtt{)} \mathtt{*} \fTHREEpx.
\end{equation}
Terms have been expanded in the table given above to clarify
how $\mathtt{ans}'$ builds up a value that is equivalent---from
the standpoint of evaluation in real arithmetic---to
$\mathtt{Prod}'\mathtt{(x)} = {\displaystyle \sum_{i=1}^k } \fipx * {\displaystyle \prod_{j \neq i} } \fjx$.
However, \eqref{ProdTraceExample} shows what happens operationally:
for a summand of the form
\[
  \fONEx * \ldots * \fimmx * \fip * \fippx * \ldots * \fkx,
\]
the product $\fONEx * \ldots * \fimmx$ is built up in variable $\mathtt{ans}$
during the first $i-1$ iterations;
then, in iteration $i$, it is multiplied by $\fip$ and added to
$\mathtt{ans}'$;
thereafter, on subsequent iterations, it resides in $\mathtt{ans}'$ 
as $\mathtt{ans}'$ is multiplied by $\fippx * \ldots * \fkx$.
\end{Exa}

For the automatic-differentiation approach, we did not really need
to make the assumption that we were given programs $\fip$ for the
functions $f_i'$, $1 \leq i \leq k$;
instead, the programs $\fip$ can be generated from the programs $\mathtt{f}_\mathtt{i}$
by applying the same statement-doubling transformation that was
applied to {\tt Prod\/}.

In languages that support operator overloading, such as C++, Ada, and
Pascal-XSC, automatic differentiation can be carried out by
defining a new data type that has fields for both the value and the
derivative, and overloading the arithmetic operators to carry out
appropriate manipulations of both fields \cite{NASC:Rall83,TOMS:Rall84},
along the lines of the definition of the C++ class {\tt FloatD},
shown in \figref{FloatD}.
\begin{figure*}
\begin{center}
{\small
  \begin{tabular}{|l|}
    \cline{1-1}
      {\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          enum ArgDesc \{ CONST, VAR \}; \\
          c\=lass FloatD \{ \+ \\
           p\=ublic: \+ \\
            float val$'$; \\
            float val; \\
            FloatD(ArgDesc,float); \-\- \\
          \}; \\
          // Constructor to convert a constant \\
          // or a value for the independent \\
          // variable to a FloatD \\
          Fl\=oatD::FloatD(ArgDesc a, float v)\{ \+ \\
            sw\=itch (a) \{ \+ \\
              ca\=se CONST: \+ \\
                val$'$ = 0.0; \\
                val = v; \- \\
              break; \\
              ca\=se VAR: \+ \\
                val$'$ = 1.0; \\
                val = v; \- \\
              break; \- \\
            \} \- \\
          \}
        \end{tabbing}
        \end{minipage}
      } \\
      {\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          Flo\=atD operator+(FloatD a, FloatD b)\{ \+ \\
            FloatD ans; \\
            ans.val$'$ = a.val$'$ + b.val$'$; \\
            ans.val = a.val + b.val; \\
            return ans; \- \\
          \}
        \end{tabbing}
        \end{minipage}
      } \\
      {\tt
        \begin{minipage}{2.0in}
        \begin{tabbing}
          Fl\=oatD operator*(FloatD a, FloatD b)\{ \+ \\
            FloatD ans; \\
            ans.val$'$ = a.val * b.val$'$ + a.val$'$ * b.val; \\
            ans.val = a.val * b.val; \\
            return ans; \- \\
          \}
        \end{tabbing}
        \end{minipage}
      } \\
    \cline{1-1}
  \end{tabular}
}
\end{center}
\caption{A differentiation-arithmetic class.}
\label{Fi:FloatD}
\end{figure*}
A class such as {\tt FloatD\/} is called a
{\em differentiation arithmetic\/} \cite{MathMag:Rall86,CASVNM:Rall90,SIAM:Rall92}.

The transformation then amounts to changing the types of each
procedure's formal parameters, local variables, and
return value (including those of the $\mathtt{f}_\mathtt{i}$).\footnote{We have
  referred to automatic differentiation as a ``program transformation,''
  which may conjure up the image of a tool that perform source-to-source rewriting fully
  automatically.  Although source-to-source rewriting is one possible embodiment,
  here the term ``transformation'' will also include the use of C++
  classes in which the arithmetic operators have been overloaded.  With
  the latter approach, rewriting might be carried out by a preprocessor,
  but might also be performed by hand, because usually only light
  rewriting of the program source text is required.
}

\begin{Exa}\label{Exa:ProdOverloadExample}
Using class {\tt FloatD\/}, the {\tt Prod\/} program of Example~\ref{Exa:ProdExample}
can be handled as follows:
\begin{center}
%\resizebox{\textwidth}{!}{
  \begin{tabular}{|@{\hspace{.3ex}}l@{\hspace{.7ex}}c@{\hspace{.8ex}}l@{\hspace{0ex}}|}
    \cline{1-3}
        {\tt float $\mathtt{f}_\mathtt{1}$(float x)\{...\}\/} & $\Rightarrow$ & {\tt FloatD $\mathtt{f}_\mathtt{1}$(const FloatD \&x)\{...\}\/} \\
        \hspace{5em}\vdots & & \hspace{5em}\vdots \\
        {\tt float $\mathtt{f}_\mathtt{k}$(float x)\{...\}\/} & $\Rightarrow$ & {\tt FloatD $\mathtt{f}_\mathtt{k}$(const FloatD \&x)\{...\}\/} \\
        {\tt
          \begin{minipage}{2.0in}
          \begin{tabbing}
            fl\=oat Prod(float x)\{ \+ \\
              float ans = 1.0; \\
              fo\=r (int i = 1; i <= k; i++)\{ \+ \\
                ans = ans * $\fix$; \- \\
              \} \\
              return ans; \- \\
            \}
          \end{tabbing}
          \end{minipage}
        }
      & $\Rightarrow$
      & {\tt
          \begin{minipage}{2.0in}
          \begin{tabbing}
            Fl\=oatD Prod(const FloatD \&x)\{ \+ \\
              FloatD ans(CONST,1.0); // ans = 1.0 \\
              fo\=r (int i = 1; i <= k; i++)\{ \+ \\
                ans = ans * $\fix$; \- \\
              \} \\
              return ans; \- \\
            \}
          \end{tabbing}
          \end{minipage}
        } \\
      & $\Rightarrow$
      & {\tt
          \begin{minipage}{2.0in}
          \begin{tabbing}
            fl\=oat Prod$'$(float x)\{ \+ \\
              FloatD xD(VAR,x); \\
              return Prod(xD).val$'$; \- \\
            \}
          \end{tabbing}
          \end{minipage}
        } \\
    \cline{1-3}
  \end{tabular}
%}
\end{center}
By changing the types of the formal parameters, local
variables, and the return values of {\tt Prod\/} and the $\mathtt{f}_\mathtt{i}$ (and making
a slight change to the initialization of {\tt ans\/} in {\tt Prod\/}), the
program now carries around derivative values (in the {\tt val\/}$'$
field) in addition to performing all of the work performed by the
original program.
Because of the C++ overload-resolution mechanism,
the $\mathtt{f}_\mathtt{i}$ procedures invoked in the fourth line of the transformed
version of {\tt Prod\/} are the {\em transformed\/} versions of the $\mathtt{f}_\mathtt{i}$
(i.e., the $\mathtt{f}_\mathtt{i}$ of type $\mathtt{FloatD} \rightarrow \mathtt{FloatD}$).

The value of {\tt Prod\/}'s derivative at {\tt v\/} is obtained by
calling $\mathtt{Prod}'\mathtt{(v)}$.
\end{Exa}

In a differentiation arithmetic, each procedure in the user's program,
such as {\tt Prod\/} and the $\mathtt{f}_\mathtt{i}$ in Example~\ref{Exa:ProdOverloadExample},
can be viewed as a box that maps two inputs to two outputs, as
depicted below:
\begin{center}
  \includegraphics[height=1.5in]{figures/compDiff-eps-converted-to.pdf}
\end{center}
In particular, in each differentiating version of a user-defined or
library procedure {\tt F\/}, the lower-right-hand output produces the
value {\tt F$'$(v)*w\/}.

An input value {\tt v\/} for the formal parameter is treated as a pair
{\tt (v,1.0)\/}.
Boxes like the one shown above ``snap together'':
when {\tt F\/} is composed with {\tt G\/} (and the input is {\tt v\/}), the
output value on the lower-right-hand side is
{\tt F$'$(G(v))*G$'$(v)\/}, which agrees with the usual
expression for the chain rule for the first-derivative operator:
\begin{center}
%\resizebox{.95\textwidth}{!}{
  \includegraphics[height=1.5in]{figures/composition-eps-converted-to.pdf}
%}
\end{center}

It should be noted that the transformation that was illustrated in
Example~\ref{Exa:ProdExample} is not fully general in that it
does not yield a procedure that can be composed with other
transformed procedures.
To create a composable transformed procedure, the transformation
would, in essence, have to make changes that mimic all of the
actions of the version created in Example~\ref{Exa:ProdOverloadExample}
using class {\tt FloatD\/}:
the procedure would have to take two arguments, {\tt x\/} and
$\mathtt{x}'$; pass these on to composable transformed versions of the
$\mathtt{f}_\mathtt{i}$; and return a pair
$\B{\mathtt{ans},\mathtt{ans}'}$, instead of $\mathtt{ans}'$ alone.

The automatic-differentiation technique summarized above
is what is known as {\em forward-mode\/} differentiation.
The availability of overloading makes it possible to implement
forward-mode automatic differentiation conveniently, by
packaging it as a differentiation-arithmetic class, as illustrated
above.
The alternative to the use of overloading is to build a
special-purpose preprocessor to carry out the statement-doubling
transformation that was illustrated in
\exrefs{ProdExample}{ProdTraceExample}.
Examples of systems that use the latter approach include ADIFOR
\cite{SC:BCGGH92,IEEE-CSE:BCKM96} and ADIC \cite{SPE:BRM97}.

When the number of independent variables is much greater than the
number of dependent variables, a different technique, {\em reverse-mode
automatic-differentiation}
\cite{BIT:Linn76,phd:Speel80,JJAM:Iri84,MPRDA:Griewank89,SIAMNews:Griewank91},
provides theoretically better performance---a greatly reduced number
of computation steps---but at the cost of the need to store
or recompute intermediate values that affect the final result
nonlinearly.
Reverse-mode automatic differentiation is explained in
\sectseqref{ComputationGraphs}{DerivativesForOtherOperators}.
In \sectref{BackpropagationForTrainingANeuralNet}, we describe how
reverse-mode automatic differentiation is the workhorse of the
\emph{backpropagation algorithm} used to train a neural network.

Forward mode and reverse mode are the endpoints of a spectrum
of algorithmic techniques;
in practice, automatic-differentiation tools
optimize runtime and memory requirements by exploiting
associativity properties of the chain rule to permit
forward mode and reverse mode to be used in different parts
of the computation \cite{SIAM:BH96}.

\section{Computation Graphs}
\label{Se:ComputationGraphs}

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{6.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=6.0cm]{figures/expression-graph-eps-converted-to.pdf}
  &
  \includegraphics[height=6.0cm]{figures/expression-graph-evaluated-eps-converted-to.pdf}
  \\
  (a) & (b)
\end{tabular}
\caption{(a) A symbolic computation graph for the evaluation of the expression $(x+4y) * ((x+y) * (y))$.
         (b) The computation graph for $(x+4y) * ((x+y) * (y))$ when $x = 3$ and $y = 5$.
  (It is recommended to view the figures in this document on-line so that the
  colors in them can be seen.)
}
\label{Fi:ExpressionGraph}
\end{figure}

In previous lectures, we described how interprocedural
dataflow-analysis problems can be formalized as path problems on graphs.
Here, we will reformulate automatic differentiation---both forward-mode
and reverse-mode---as a path problem.
To do so, we introduce the notion of a \textit{computation graph}
\cite{Backprop}.
A computation graph is a labeled directed-acyclic graph (DAG).
We start with the---in general---multi-source/multi-sink DAG that
consists of all operations performed during execution, connected
by edges that represent how outputs from one operation were used
as inputs to another operation.
\figref{ExpressionGraph}(a) shows an example of a computation graph
for the following function:
\begin{equation}
  \label{Eq:ExampleFunction}
  h(x,y) = (x+4y) * ((x+y) * (y)).
\end{equation}
We will use \eqref{ExampleFunction} as a running example in the rest
of this document.
Note that by expanding the right-hand-side expression,
$h(x,y)$ equals $x^2y + 5xy^2 + 4y^3$ (i.e., it is
a multivariate polynomial).

To formulate automatic differentiation as a path problem, we have
to define the values that label a computation graph's edges
and to define two operators, \emph{extend} and \emph{combine}---as
we did for shortest-distance problems and dataflow-analysis problems.
The path problem we need to define here is relatively straightforward:
the extend operator is \emph{multiplication} and the combine operator is
\emph{addition}.
Note that unlike in the dataflow-analysis case, the combine operator is
not idempotent:
for instance, $1+1 = 2$ and not $1$!
However, the absense of idempotence does not create difficulties
for automatic differentiation because we only need to solve
path problems on a \emph{DAG};
they can always be solved using a traversal of the DAG in topological order.

The edges are labeled according to the following scheme:

\noindent
\begin{center}
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{10.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=1.2in]{figures/annotation-plus-eps-converted-to.pdf}
  &
  \includegraphics[height=1.2in]{figures/annotation-times-eps-converted-to.pdf}
\end{tabular}
\end{center}

\noindent
The labeling schema given above is really a \emph{symbolic} version
of the method used to construct a computation graph;
\figref{ExpressionGraph}(a) is an example of a symbolic
computation graph.
In an actual computation graph, we will have specific values for
$x$ and $y$, as well as for each of the subexpressions used as labels:
$(x+y)y$, $x + 4y$, and $x+y$.
\figref{ExpressionGraph}(b) shows the computation graph for
$(x+4y) * ((x+y) * (y))$ when $x = 3$ and $y = 5$.

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{6.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=6.0cm]{figures/partial-x-forward-eps-converted-to.pdf}
  &
  \includegraphics[height=6.5cm]{figures/partial-x-reverse-eps-converted-to.pdf}
  \\
  (a) & (b)
  \\
  \\
  \includegraphics[height=6.0cm]{figures/partial-y-forward-eps-converted-to.pdf}
  &
  \includegraphics[height=6.0cm]{figures/partial-y-reverse-eps-converted-to.pdf}
  \\
  (c) & (d)
\end{tabular}
\caption{The computation of partial derivatives of $h(x,y)$,
  when $x = 3$ and $y = 5$, via forward and backward single-source/single-target
  path problems.
  (a) $\frac{\partial((x+4y) * ((x+y) * (y)))}{\partial x} = 155$ via a forward path problem;
  (b) $\frac{\partial((x+4y) * ((x+y) * (y)))}{\partial x} = 155$ via a backward path problem;
  (c) $\frac{\partial((x+4y) * ((x+y) * (y)))}{\partial y} = 459$ via a forward path problem;
  (d) $\frac{\partial((x+4y) * ((x+y) * (y)))}{\partial y} = 459$ via a backward path problem;
}
\label{Fi:PartialDerivatives}
\end{figure}

\section{Computing Partial Derivatives Using a Computation Graph}
\label{Se:ComputingPartialDerivativesUsingAComputationGraph}

The reason why computation graphs are interesting is because it is
possible to compute partial derivatives by solving
single-source/single-target path problems on a computation graph.
\figref{PartialDerivatives} returns to the example of
$h(x,y) = (x+4y) * ((x+y) * (y))$, and illustrates how
$\frac{\partial h(x,y)}{\partial x}$ and $\frac{\partial h(x,y)}{\partial y}$,
when $x = 3$ and $y = 5$, can be computed via forward and backward
single-source/single-target path problems.

The calculations for $\frac{\partial h}{\partial x}$ depicted in
\figref{PartialDerivatives}(a) and \figref{PartialDerivatives}(b)
involve summing over all paths from $x$ to the root of the
computation graph for $h$.
This sum can be computed in two directions (and the same result is obtained
via either method):
in \textit{forward-mode automatic differentiation}, the path-value
computations proceed from the leaf $x$ to the root of the computation graph.
(Forward-mode differentiation tracks how the value at each node
is affected by a change in the value of the input $x$.)
In contrast, with \textit{reverse-mode automatic differentiation} the
path-value computations proceed from the root of the computation graph
to the leaf $x$.
(Reverse-mode differentiation tracks how the output value of the
function is affected by a change in the value at a given node in the
computation graph.)

For example, the green values shown in
\figref{PartialDerivatives}(a) and \figref{PartialDerivatives}(b)
indicate the intermediate path-values computed for nodes that are
relevant to the final calculation during the forward and backward
traversals, respectively.
Because we are working with a DAG, the values can be computed
in linear time in the size of the relevant portion of the computation
graph by traversing that portion of the computation graph in
topological order.

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{10.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=2.4in]{figures/inductive-case-plus-eps-converted-to.pdf}
  &
  \includegraphics[height=2.4in]{figures/inductive-case-times-eps-converted-to.pdf}
  \\
  (a) & (b)
\end{tabular}
\caption{The essence of the induction step in an inductive proof of
  correctness of the method for computing a partial derivative by
  solving a forward path problem (forward-mode AD).
  Diagrams (a) and (b) show how information is propagated from
  children to their parent in the case of a $+$-node and a $\times$-node,
  respectively.
}
\label{Fi:InductiveStep}
\end{figure}

In the case of forward mode AD, the diagrams in \figref{InductiveStep} depict the
main elements of the induction step in an inductive proof of
correctness of the (forward) path-problem method.
The induction hypothesis is that the single-source/single-target
forward path problem in each of the shaded regions computes the
indicated partial derivatives.
The continuation of the path computations in a $(+,\times)$ path problem
yield the respective indicated values
$1 * \frac{\partial f}{\partial x} + 1 * \frac{\partial g}{\partial x}$ and
$g * \frac{\partial f}{\partial x} + f * \frac{\partial g}{\partial x}$.

However, if we know that the forward-path-problem method is a correct
method for computing partial derivatives, then the correctness of the
backward-path-problem method follows immediately from the associativity
and commutativity of $+$ and $\times$.

\begin{Exa}\label{Exa:PartialDerivatives}
Let us now return to \eqref{ExampleFunction} and show that the values
computed in \figref{PartialDerivatives}(a) and \figref{PartialDerivatives}(b)
are correct by computing the partial derivative by the more usual rules:
\begin{align*}
  \left.\frac{\partial h(y,x)}{\partial x}\right|_{x=3,y=5}
       & = \left.\frac{\partial (x^2y + 5xy^2 + 4y^3)}{\partial x}\right|_{x=3,y=5} \\
       & = \left.2xy + 5y^2\right|_{x=3,y=5} \\
       & = 30 + 125 \\
       & = 155
\end{align*}

Another way to check our work is to solve the path problems from
\figref{PartialDerivatives}(a) and \figref{PartialDerivatives}(b)
symbolically, using \figref{ExpressionGraph}(a), and show that
both methods produce the answer $2xy + 5y^2$.
\begin{description}
  \item [{\it Forward path problem:}]
    \begin{align*}
      \frac{\partial h(y,x)}{\partial x}
          & = 1 * (x + y)y + 1 * y * (x + 4y) \\
          & = xy + y^2 + yx + 4y^2) \\
          & = 2xy + 5y^2
    \end{align*}
  \item [{\it Backward path problem:}]
    \begin{align*}
      \frac{\partial h(y,x)}{\partial x}
          & = (x + y)y * 1 + (x + 4y) * y * 1 \\
          & = xy + y^2 + xy + 4y^2 \\
          & = 2xy + 5y^2
    \end{align*}
\end{description}
\end{Exa}

\begin{figure}
\centering
  \includegraphics[height=9.0cm]{figures/forward-vs-reverse-eps-converted-to.pdf}
\caption{Depiction of the cost difference between forward-mode and
  reverse-mode automatic differentiation, depending on what set of
  partial derivatives is desired.
}
\label{Fi:ForwardVsReverse}
\end{figure}

\subsection*{Forward-Mode Versus Reverse-Mode}

At first glance, one may think because forward-mode and reverse-mode
automatic differentiation compute the same answer, neither method
is to be preferred over the other.
However, the two methods may have substantially different costs, depending
on the number of inputs and outputs of the function.

To understand the cost difference between forward-mode and
reverse-mode automatic differentiation, consider a multivariate
function $\vec{h}(\vec{x})$ that maps a vector of $n$ inputs
$\B{x_1, x_2, \ldots, x_n}$ to a vector of $m$ outputs
$\B{h_1(\vec{x}), h_2(\vec{x}), \ldots, h_m(\vec{x})}$, as depicted
in \figref{ForwardVsReverse}.
\figref{ForwardVsReverse} illustrates how forward-mode and reverse-mode
can have much different costs, depending on what set of
partial derivatives is desired:
\begin{itemize}
  \item
    One extreme is when we want to compute the set
    $\left\{ \frac{\partial h_k}{\partial x_j} \mid 1 \leq k \leq m \right\}$
    of \emph{all} partial derivatives of $\{ h_k(\vec{x}) \mid 1 \leq k \leq m \}$
    with respect to \emph{some} input $x_j$.
    In effect, $n = 1$ (although we must have available the results of
    any sub-computation that can affect the value of any of the $h_k(\vec{x})$).
    In this case, the set of partial derivatives should be computed by solving a
    single-source/multi-target path problem, which can be performed
    using a single forward pass over the computation graph.
    This computation corresponds to forward-mode automatic differentiation,
    and is depicted in blue in \figref{ForwardVsReverse}.
  \item
    Another extreme is when we want to compute the set
    $\left\{ \frac{\partial h_i}{\partial x_k} \mid 1 \leq k \leq n \right\}$
    of partial derivatives of \emph{some} $h_i(\vec{x})$ with
    respect to \emph{all} inputs $\{ x_k \mid 1 \leq k \leq n \}$.
    In effect, $m = 1$ (and we really only need the computation graph
    for $h_i(\vec{x})$).
    In this case, the set of partial derivatives should be computed by
    solving a multi-source/single-target path problem, which can be
    performed using a single backward pass over the computation
    graph.
    This computation corresponds to reverse-mode automatic
    differentiation, and is depicted in purple in
    \figref{ForwardVsReverse}.
\end{itemize}
More generally, we might want the set
$\left\{ \frac{\partial h_k}{\partial x_l} \mid k \in K \subseteq \{1, \ldots, m\}
                                          \land l \in L \subseteq \{1, \ldots, n\} \right\}$.
There is a trade-off between performing $|L|$ forward passes or
$|K|$ backward passes.
Each forward pass computes $\left\{ \frac{\partial h_k}{\partial x_j} \mid k \in K \right\}$
for some $j \in L$.
Each backward pass computes $\left\{ \frac{\partial h_i}{\partial x_l} \mid l \in L \right\}$
for some $i \in K$.
In this case, reverse-mode automatic differentiation is preferred
when $|L| \gg |K|$.
Forward-mode automatic differentiation is preferred when $|L| \ll |K|$.

As we will see in \sectref{BackpropagationForTrainingANeuralNet},
when training a neural network, the number of inputs is typically
large, and the number of outputs is $1$.
Consequently, the backpropagation algorithm for neural-net training
employs reverse mode.

\section{Derivatives for Other Operators}
\label{Se:DerivativesForOtherOperators}

In this section, we describe how to construct the computation graph
when the computation contains operators other than addition and multiplication.
In general, for an occurrence of an $n$-ary operator, the label on the
$j^{\textit{th}}$ argument is the value of the partial derivative
$\frac{\partial \textit{op}(f_1, \ldots, f_j, \ldots, f_n)}{\partial f_j}$:

\noindent
\begin{center}
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=.8in]{figures/n-ary-op-eps-converted-to.pdf}
\end{tabular}
\end{center}

\noindent
It is easy to see that the two labeling schemes for addition and
multiplication depicted at the end of \sectref{ComputationGraphs}
fit the general pattern.
We now consider several other examples.

\paragraph*{Square Root.}

The square-root operator $\sqrt x \eqdef x^{\frac{1}{2}}$ is a unary operator;
consequently, there is only one incoming edge to a square-root node.
The value on the edge is the value of the partial derivative
$\frac{\partial \sqrt x}{\partial x}$
$= \frac{\partial (x^{\frac{1}{2}})}{\partial x}$
$= \frac{1}{2} x^{-\frac{1}{2}}$
$= \frac{1}{2\sqrt x}$.
Consequently, for an occurrence of the square-root operator
that produces the value $v$, the label that should be placed
on the incoming edge is $\frac{1}{2 v}$.

\noindent
\begin{center}
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=.8in]{figures/sqrt-clipped.pdf}
\end{tabular}
\end{center}

\noindent
Note that the square-root operator is different from addition and multiplication
in that the value on the edge is a function of the \emph{output}
of the operator.

\paragraph*{Sigmoid.}

One of the operations used in neural nets is the sigmoid operation, which
is also unary:
\[
  \sigma(x) \eqdef \frac{1}{1 + e^{-x}}.
\]
The value on the edge is the value of the partial derivative
{\small
\begin{align*}
  \frac{\partial \sigma(x)}{\partial x}
      & = \frac{\partial {(1+e^{-x})}^{-1}}{\partial x} \\
      & = -1(1+e^{-x})^{-2} e^{-x} (-1) \\
      & = \frac{e^{-x}}{(1+e^{-x})^{2}} \\ 
      & = \left(\frac{1}{1+e^{-x}}\right) \left(\frac{e^{-x}}{1+e^{-x}}\right) \\
      & = \left(\frac{1}{1+e^{-x}}\right) \left(\frac{1+e^{-x}}{1+e^{-x}} - \frac{1}{1+e^{-x}}\right) \\
      & = \sigma(x) (1 - \sigma(x))
\end{align*}}
Consequently, for an occurrence of the sigmoid operator
that produces the value $v$, the label that should be placed
on the incoming edge is $v(1 - v)$.

\noindent
\begin{center}
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=.8in]{figures/sigmoid-eps-converted-to.pdf}
\end{tabular}
\end{center}


\paragraph*{Squares of Differences.}

Consider the following 2-argument function used for computing one term in a sum of squares of
differences:
\[
  \lambda x,v. \frac{1}{2}(x-v)^2.
\]
In the context of training a neural net, the first argument will be the value computed
by one output node of the neural net.
The second argument will be the expected value of a training example, and
thus is considered to be a constant.
Consequently, we are only interested in the value on the left edge, which
is the value of the partial derivative
\[
  \frac{\partial \frac{1}{2}(x-v)^2}{\partial x} = x - v.
\]

\noindent
\begin{center}
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=1.3in]{figures/squared-difference-eps-converted-to.pdf}
\end{tabular}
\end{center}

\section{Backpropagation for Training a Neural Net}
\label{Se:BackpropagationForTrainingANeuralNet}

In this section, we explain how reverse-mode automatic differentiation
forms the basis of the backpropagation algorithm for training a neural net:
in backpropagation, the weight adjustments computed for a neural net,
a given input, and an expected output are found by performing
reverse-mode automatic differentiation on the computation graph for
(i) the neural net, along with (ii) additional computation-graph elements for the sum of
the squared differences of the neural net's output with the expected output.

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{4.0ex}}l@{\hspace{2.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=5.0cm]{figures/neural-net-eps-converted-to.pdf}
  &
  \begin{minipage}[b]{2.5cm}
    training layer \\ \\ \\
    output layer \\ \\ \\
    hidden layer \\ \\ \\
    input layer \\
  \end{minipage}
  &
  \includegraphics[height=6.9cm]{figures/neural-net-training-eps-converted-to.pdf}
  \\
  (a) & & (b)
\end{tabular}
\caption{(a) A neural net with one hidden layer.
         (b) Circuit used for training a neural net.
}
\label{Fi:NeuralNet}
\end{figure}

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{4.0ex}}l@{\hspace{2.0ex}}c@{\hspace{0ex}}}
  \includegraphics[height=4.9cm]{figures/neuron-eps-converted-to.pdf}
\end{tabular}
\caption{Structure of a node in a neural network.}
\label{Fi:Neuron}
\end{figure}

A neural network has an input layer, an output layer, and at least one
hidden layer in the middle, as shown in \figref{NeuralNet}(a).
When one looks into one of the nodes in the network, one typically
has what is depicted in \figref{Neuron}:
\begin{itemize}
  \item
    a node has some number of inputs $i_1, \ldots, i_k$, which are multiplied
    respectively by the weights $w_1, \ldots, w_k$
  \item
    these products are then summed
  \item
    a sigmoid function is applied to the sum
  \item
    the output of the sigmoid function goes off to the next layer
    in the network---in general, to multiple nodes in the next layer
\end{itemize}

A neural network is trained by a ``hill-climbing'' procedure.
One repeatedly presents it with training examples---i.e., pairs of
vectors $\B{\vec{x},\vec{v}}$, where $\vec{v}$ represents the desired
output for input $\vec{x}$---drawn from some suite of training data.
For understanding how training works, it is useful to think of the
neural net as being equipped with an additional layer that computes a
cost $\Cost(\vec{y},\vec{v})$ by computing the sum of the squared
differences of the network's actual output for input $\vec{x}$ (i.e.,
$\vec{y}$) and the expected output $\vec{v}$---see
\figref{NeuralNet}(b).
The use of the sum-of-squared differences
\[
  \sum_{1 \le j \le m} \frac{1}{2}(y_j - v_j)^2
\]
helps to accommodate noise in the training data.
(The training data might contain mistakes due to human error, or it might
consist of data observed from the environment via sensors, which will introduce
measurement error.)

For each $\B{\vec{x},\vec{v}}$ pair, the weights in each of the nodes are adjusted by
\begin{enumerate}
  \item
    using reverse-mode automatic differentiation to compute $\frac{\partial \Cost}{\partial w_i}$,
    for each $w_i$ that is a weight component of one of the hidden or output nodes
    in the neural network
  \item
    making the adjustment $w_i := w_i + \eta \frac{\partial \Cost}{\partial w_i}$,
    where $\eta$ is a small constant that controls the rate of hill-climbing.
\end{enumerate}
Training continues until the weight adjustments fall below some desired threshold.

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[width=.6\textwidth]{figures/neural-net-weights-as-inputs-clipped.pdf}
\end{tabular}
\caption{A neural network has three kinds of inputs: (i) the input
  data $\vec{x}$, (ii) the parameters $\vec{w}$, and (iii) the desired
  output $\vec{v}$.
  To perform back-propagation, one needs to compute the partial derivatives
  $\left\{ \frac{\partial \Cost}{\partial w_i} \right\}$.
  (Adapted from a slide in a POPL 2022 talk by Faustyna Krawiec.)
  }
\label{Fi:AbstractExample}
\end{figure}

Note that when training a neural network using automatic differentation, it is the set
of weights $\{ w_i \}$ that are considered to be the inputs of the computation graph
with respect to which we want to compute the partial derivatives
$\left\{ \frac{\partial \Cost}{\partial w_i} \right\}$ (see \figref{AbstractExample}).
The values from $\vec{x}$ are treated as constants for the purposes of a given round
of training.
Moreover, $\Cost$ is the only output, and thus the number of outputs is $1$.
Consequently, automatic differentiation is performed using reverse mode.

\begin{figure}
\centering
\begin{tabular}{@{\hspace{0ex}}c@{\hspace{0ex}}}
  \includegraphics[width=.6\textwidth]{figures/full-example-eps-converted-to.pdf}
\end{tabular}
\caption{A back-propagation example in schematic form.}
\label{Fi:FullExample}
\end{figure}

\figref{FullExample} shows a full neural-net-training example in schematic form.
The figure shows the computation graph for a two-layer
neural net, including the training layer.
The edges of the graph are annotated with the edge-weights used in the path problem
for automatic differentiation.
With reverse-mode automatic differentiation, the desired set of partial derivatives
$\left\{ \frac{\partial \Cost}{\partial w_i} \right\}$ can be computed using a single backward
pass over the computation graph.

It is instructive to consider in more detail a couple of the partial derivatives
that are computed.
For instance, there is only one path from $\Cost$ to $w_4$,
and thus $\frac{\partial \Cost}{\partial w_4}$ can be computed
as follows:
\[
  \begin{array}{@{\hspace{0ex}}r@{\hspace{1.0ex}}c@{\hspace{1.0ex}}l@{\hspace{0ex}}}
    \frac{\partial \Cost}{\partial w_4} & = & 1 \times (y_2 - v_2) \times y_2(1-y_2) \times 1 \times h_1 \\
                                        & = & (y_2 - v_2) y_2(1-y_2) h_1
  \end{array}
\]
In contrast, there are multiple paths from $\Cost$ to $w_{14}$:
\[
  \begin{array}{@{\hspace{0ex}}r@{\hspace{1.0ex}}c@{\hspace{1.0ex}}l@{\hspace{0ex}}}
    \frac{\partial \Cost}{\partial w_{14}}
    & = &
      \left(\begin{array}{@{\hspace{0ex}}c@{\hspace{1.0ex}}l@{\hspace{0ex}}}
          & 1 \times (y_1 - v_1) \times y_1(1-y_1) \times 1 \times w_3 \times h_3(1-h_3) \times 1 \times x_1 \\
        + & 1 \times (y_2 - v_2) \times y_2(1-y_2) \times 1 \times w_6 \times h_3(1-h_3) \times 1 \times x_1 \\
        + & 1 \times (y_3 - v_3) \times y_3(1-y_3) \times 1 \times w_9 \times h_3(1-h_3) \times 1 \times x_1
      \end{array}\right) \\
    & = &
      \left(\begin{array}{@{\hspace{0ex}}c@{\hspace{1.0ex}}l@{\hspace{0ex}}}
          & (y_1 - v_1)y_1(1-y_1)w_3 \\
        + & (y_2 - v_2)y_2(1-y_2)w_6 \\
        + & (y_3 - v_3)y_3(1-y_3)w_9
      \end{array}\right) h_3(1-h_3)x_1
  \end{array}
\]
These results are then used to adjust the values of $w_4$ and $w_{14}$ by making the adjustments
\begin{align*}
  w_4    &= w_4 + \eta \frac{\partial \Cost}{\partial w_4} \\
  w_{14} &= w_{14} + \eta \frac{\partial \Cost}{\partial w_{14}}
\end{align*}

{\small
  \bibliography{cs701,df,logic,diff}
  \bibliographystyle{plain}
}

\end{document}
