\documentclass[12pt]{article}
\usepackage{epsfig}
\usepackage{url}
\usepackage{verbatim}

\newcommand{\yy}[1]{\textit{todo: #1}}
\newcommand{\ct}[1]{\textit{cite: #1}}
\newcommand{\cd}[1]{\texttt{#1}}
\newcommand{\trm}[1]{\textit{#1}}
\newcommand{\vek}[1]{\textbf{#1}}

\newcommand{\loss}{\textit{L}}
\newcommand{\deriv}[2]{{\frac{\displaystyle d{}#1}{\displaystyle d{}#2}}}

\title{Automatic Reverse-Mode Differentiation}

\author{William Cohen\\ Modified by xxx}
\date{Out xxx\\Due xxx via Blackboard}

\begin{document}

\maketitle

\section{Background: Automatic Differentiation}

\subsection{Why automatic differentiation is important}

Most neural network packages (e.g., Torch or Theano) don't require a
user to actually derive the gradient updates for a neural model.
Instead they allow the user to define a \trm{model} in a ``little
language'', which supports common neural-network operations like
matrix multiplication, ``soft max'', etc, and automatically derive the
gradients.  Typically, the user will define a \trm{loss function}
$\loss$ in this sublanguage: the loss function takes as inputs the
training data, current parameter values, and hyperparameters, and
outputs a scalar value. From the definition of $\loss$, the system
will compute the partial derivative of the loss function with respect
to every parameter, using these gradients, it's straightforward bto
implement gradient-based learning methods.

Going from the definition of $\loss$ to its partial derivatives is
called \trm{automatic differentiation}.  In this assignment you will
start with a simple automatic differentiation system written in
Python, and use it to implement a neural-network package.

\subsection{Wengert lists}

Automatic differentiation proceeds in two stages.  First, function
definitions $f(x_1,\ldots,x_k)$ in the sublanguage are converted to a
format called a \trm{Wengert list}.  The second stage is to evaluate
the function and its gradients using the Wengert list.

A \trm{Wengert list} defines a function $f(x_1,\ldots,x_k)$ of some
inputs $x_1,\ldots,x_k$.  Syntactically, it is just a list of
assignment statements, where the right-hand-side (RHS) of each
statement is very simple: a call to a function $g$ (where $g$ is one
of a set of primitive functions $G=\{g_1,\ldots,g_\ell\}$ supported by
the sublanguage), where the arguments to $g$ either inputs
$x_1,\ldots,x_k$, or the left-hand-side (LHS) of a previous
assignment.  (The functions in $G$ will be called \trm{operators}
here.)  The output of the function is the LHS of the last item in the
list.  For example, a Wengert list for
\begin{equation} \label{eq:2x1Plusx2}
f(x_1,x_2) \equiv (2 x_1 + x_2)^2
\end{equation} 
with the functions $G = \{ \cd{add},\cd{multiply},\cd{square} \}$ might be
\begin{eqnarray*}
z_1 & = & \cd{add}(x_1,x_1) \\
z_2 & = & \cd{add}(z_1,x_2) \\
f & = & \cd{square}(z_2)
\end{eqnarray*}
A Wengert list for
\[ f(x) \equiv x^3
\] 
might be
\begin{eqnarray*}
z_1 & = & \cd{multiply}(x,x) \\
f & = & \cd{multiply}(z_1,x)
\end{eqnarray*}

The set of functions $G$ defines the sublanguage.  It's convenient if
they have a small fixed number of arguments, and are differentiable
with respect to all their arguments.  But there's no reason that they
have to be scalar functions!  For instance, if $G$ contains the
appropriate matrix and vector operations, the loss for logistic
regression (for example $\vek{x}$ with label $y$ and weight matrix
$W$) could be written as
\[ f(\vek{x},\vek{y},W) \equiv 
      \cd{crossEntropy}( \cd{softmax}(\vek{x} \cdot W ), \vek{y}) + \cd{frobeniusNorm}(W)
\] 
and be compiled to the Wengert list
\begin{eqnarray*}
 \vek{z}_1 & = & \cd{dot}(\vek{x},W) \\
 \vek{z}_2 & = & \cd{softmax}(\vek{z}_1) \\
 z_3 & = & \cd{crossEntropy}(\vek{z}_2,\vek{y}) \\
 z_4 & = & \cd{frobeniusNorm}(W) \\
 f & = & \cd{add}(z_3,z_4)
\end{eqnarray*}
This implements $k$-class logistic regression if $\vek{x}$ is a
$d$-dimensional row vector, \cd{dot} is matrix product, $W$ is a $d
\times k$ weight matrix, and 
\begin{eqnarray*}
 \cd{softmax}(\langle a_1,\ldots,a_d\rangle) & \equiv &
   \langle \frac{e^{a_1}}{\sum_{i=1}^d e^a_i}, \ldots, \frac{e^{a_d}}{\sum_{i=1}^d e^a_i} \rangle  \\
 \cd{crossEntropy}(\langle a_1,\ldots,a_d\rangle,\langle b_1,\ldots,b_d\rangle) & \equiv &
    - \sum_{i=1}^d a_i \log b_i \\
 \cd{frobeniusNorm}(A) & \equiv & \sqrt{ \sum_{i=1}^d \sum_{j=1}^k a_{i,j}^2 }
\end{eqnarray*}

\subsection{Using a Wengert list}

We'll first discuss how to use a Wengert list, and then below, discuss
how to construct one.

Given a Wengert list for $f$, it's obvious how to evaluate $f$: just
step through the assignment statements in order, computing each value
as needed.  To be perfectly clear about this, the procedure is as
follows.  We will encode each assignment in Python as a nested tuple
\[ (z, g, (y_1,\ldots,y_k))
\]
where $z$ is a string that names the LHS variable, $g$ is a string
that names the operator, and the $y_i$'s are strings that name the
arguments to $g$.  So the list for the function of
Equation~\ref{eq:2x1Plusx2} would be encoded in python as
\begin{verbatim}
[ ("z1", "add", ("x1","x2")), 
  ("z2", "add", ("z1","x2")), 
  ("f", "square", ("z_2")) ]
\end{verbatim}
We also store functions for each operator in a Python dictionary
\cd{G}:
\begin{verbatim}
G = { "add" : lambda a,b: a+b, 
       "square": lambda a:a*a }
\end{verbatim}
Before we evaluate the function, we will store the parameters in a
dictionary \cd{val}: e.g., to evalute $f$ and $x_1=3$, $x_2=7$
we will initialize \cd{val} to 
\begin{verbatim}
val = { "x1" : 3, "x2" : 7 }
\end{verbatim}

The pseudo-code to evaluate $f$ is:

\begin{tabbing}1234\=1234\=\kill
\cd{def eval($f$)} \\
\> initialize \cd{val} to the inputs at which $f$ should be evaluated\\
\> \cd{for} $(z,g,(y_1,\ldots,y_k))$ in the list:\\
\>     \> \cd{op = G[g]} \\
\>     \> \cd{val[z] = op(val[$y_1$],\ldots,val[$y_k$])}\\
\> \cd{return} the last entry stored in \cd{val}.
\end{tabbing}

\textit{Some Python hints:} (1) to convert $(y_1,\ldots,y_k)$ to
\cd{(val[$y_1$],\ldots,val[$y_k$])} you might use Python's \cd{map}
function. (2) If \cd{args} is a length-2 Python list and \cd{g} is a
function that takes two arguments (like \cd{G[``add'']} above) then
\cd{g(*args)} will call \cd{g} with the elements of that list as the
two arguments to \cd{g}.

To differentiate, we will use a generalization of backpropagation
(backprop). We'll assume that \cd{eval} has already been run and
\cd{val} has been populated, and we will compute, in \emph{reverse}
order, a value \cd{delta}($z_i$) for each variable $z_i$ that appears
in the list.  We initialize \cd{delta} by setting $\cd{delta}(f)=1$,
(where $f$ is the string that names the function output).  

Informally you can think of \cd{delta}($z_i$) as the ``sensitivity''
of $f$ to the variable $z_i$, , at the point we're evaluating $f$
(i.e., the point $a$ that corresponds to the initial dictionary
entries we stored in \cd{val}). Here $z_i$ can be intermediate
variable or an input.  If it's an input $x$ that we're treating as a
parameter, \cd{delta} is the gradient of the cost function, evaluated
at $a$: i.e., $\cd{delta}(x) = \deriv{f}{x}(a)$.  

To compute these sensitivities we need to ``backpropagate'' through
the list: when we encounter the assignment \((z,g,(y_1,\ldots,y_k))\)
we will use \cd{delta}($z$) and the derivatives of $g$ with respect to
its inputs to compute the sensitivities of the $y$'s.  We will store
derivatives for each operator in a Python dictionary \cd{DG}. 

Note that if $g$ has $k$ inputs, then we need $k$ partial derivatives,
one for each input, so the entries in \cd{DG} are \emph{lists} of
functions. For the functions used in this example, we'll need
these entries in DG.
\begin{verbatim}
DG = { "add" : [ (lambda a,b: 1),  (lambda a,b: 1) ],  
       "square": [ lambda a:2*a ] }
\end{verbatim}
To figure out these functions we used some high school calculus rules:
$\deriv{}{x}(x+y)=1$, $\deriv{}{y}(x+y)=1$ and $\deriv{}{x}(x^2)=2x$.

Finally, the pseudo-code to compute the deltas is below.  Note that we
don't just store values in \cd{delta}: we accumulate them additively.

\begin{tabbing}1234\=1234\=1234\=1234\=\kill
\cd{def backprop($f$,val)} \\
\> initialize \cd{delta}: \cd{delta}[$f$] = 1\\
\> \cd{for} \cd{(z,g,($y_1,\ldots,y_k$))} in the list, in reverse order:\\
\>     \> \cd{for} $i=1,\ldots,k$:\\
\>     \>    \> \cd{op$_i$ = DG[g][$i$]} \\
\>     \>    \> \cd{if delta[$y_i$]} is not defined set \cd{delta[$y_i$] = 0} \\
\>     \>    \> \cd{delta}[$y_i$] = \cd{delta}[$y_i$] + \cd{delta[z] * op$_i$(val[$y_1$],\ldots,val[$y_k$])}\\
\end{tabbing}

What's going on here? Let's simplify for a minute and assume that the
list is of the form
\begin{eqnarray*}
 z_1 & = &  f_1(z_0) \\
 z_2 & = &  f_2(z_1) \\
 & \ldots & \\
 z_m & = &  f_m(z_{m-1})
\end{eqnarray*}
so $f = z_m = f_m(f_{m-1}(\ldots f_1(z_0) \ldots))$.  We'll assume we
can compute the $f_i$ functions and their derivations $f'_i$.  We know
that one way to find $\deriv{z_m}{z_0}$ would be to repeatedly use the
chain rule:
\begin{eqnarray*}
\deriv{z_m}{z_0} & = & \deriv{z_m}{z_{m-1}} \deriv{z_{m-1}}{z_0} \\
                 & = & \deriv{z_m}{z_{m-1}} \deriv{z_{m-1}}{z_{m-2}} \deriv{z_{m-2}}{z_{0}} \\
 \ldots \\
                 & = &  \deriv{z_m}{z_{m-1}} \deriv{z_{m-1}}{z_{m-2}} \ldots \deriv{z_{1}}{z_{0}}
\end{eqnarray*}
Let's take some time to unpack what this means.  To simplify notation,
let $h_{i,j}$ be the function $\deriv{z_i}{z_j}$.  Notice that there
are $m^2$ of these $h_{i,j}$ functions---quite a few!  We'll look at
evaluating $h_{m,0}$ as some point $z_0=a$ (say $z_0=53.4$). Again to
simplify, define 
\begin{eqnarray*}
a_1 & = & f_1(a) \\
a_2 & = & f_2(f_1(a)) \\
 & \ldots & \\
a_m & = & f_{m}(f_{m-1}(f_{m-2}(\ldots f_1(a) \ldots)))
\end{eqnarray*}
When we write
\[
\deriv{z_m}{z_0}  =  \deriv{z_m}{z_{m-1}} \deriv{z_{m-1}}{z_0}
\]
We mean that for all $a$
\[
h_{m,0}(a) = f'_{m}(a_m) \cdot h_{m-1,1}(a)
\]
Likewise
\begin{eqnarray*}
\deriv{z_m}{z_0} & = & \deriv{z_m}{z_{m-1}} \deriv{z_{m-1}}{z_{m-2}} \ldots \deriv{z_{1}}{z_{0}}
\end{eqnarray*}
means that 
\[
h_{m,0}(a) = f'_{m}(a_m) \cdot f'_{m-1}(a_{m-1})\cdots f'_{2}(a_1) \cdot f'_{1}(a)
\]
which conveniently can be coded up to actually evaluate $h_{m,0}$
given procedures $f'_i$ to compute the local function derivatives.

\yy{mapping to code, $z_i$ is a string, $a_i$ is val[$z_i$], and
  $f'_{m}(a_m) \cdots f'_{i}(a_{i})$ is delta[$z_i$]}

Let's see how this works for Equation~\ref{eq:2x1Plusx2}. 


\subsection{Constructing Wengert lists}

\input{example}


\section{Assignment}

\yy{more hints: return delta[z]*d, to you can do the matrix-vector
  thing; sometimes output is used in derivatives.}

\section{Data}

\section{Deliverables}

\section{Marking breakdown}

\end{document}

