\documentstyle[fullpage]{article}
\begin{document}

\newcommand{\rR}{I \!\! R}
\newcommand{\AD}{{\sl Algorithmic Differentiation }}

\begin{center}
{\LARGE\bf GRADIENT}\\[0.4cm] 
{\LARGE\bf Algorithmic Differentiation in Maple} \\
\bigskip
 Michael B.~Monagan\footnote{Institut f\"ur Wissenschaftliches Rechnen,
Symbolic Computation Group, ETH Z\"urich\\CH--8092 Z\"urich, Switzerland}\\
 Walter M.~Neuenschwander
\end{center}

\begin{abstract}
Many scientific applications require computation of the derivatives of a function 
$f : \rR^n \to \rR^m$ as well as the function values of $f$ itself.
All computer algebra systems can differentiate functions represented by formulae.
But not all functions can be described by formulae.
And formulae are not always the most efficient means for representing
functions and derivatives.

\verb/GRADIENT/ is a Maple \cite{Mapl91} routine that accepts a 
Maple-procedure for the
computation of a function and produces a new procedure for the computation
of its derivatives. There exists already a large amount of such AD--Tools
written in different programming languages. \verb/GRADIENT/ is one of the
first that's completely implemented in a computer algebra system like Maple.
Based on the
\verb/GRADIENT/ function it is possible to generate Maple--procedures
for the computation of the Jacobian matrix, the Hessian matrix and
coefficients of multivariate Taylor series.\\[0.3cm]
{\bf Keywords}: Algorithmic Differentiation, Chain Rule, Computer Algebra,
Gradients, Jacobians, Hessians, Taylor Coefficients
\end{abstract}

\section{Introduction}
The techniques of \AD are increasingly used in
numerical computation and applied-analytical methods. The reason for this
isn't purely for mathematical results. The usual methods based on
{\sl Divided Differences} are very inaccurate. Especially derivatives of 
order 3 and higher are completely useless. The other fact is that \AD
is much more general than other methods. The very first approaches in \AD
were done more than 30 years ago \cite{Beda59,Weng64}. Since then various 
implementations of \AD have have been done. The paper of D.~Juedes 
\cite{Jued91} gives a summary of all tools developed up till now.\\
What do we understand by \AD \/? What we 
want to do is to compute the partial derivatives of a function $f$ with
respect to an arbitrary set of variables. The function doesn't have to
be defined by a closed formula nor by a straight-forward program.
Actually we want to differentiate functions which are represented by a 
general algorithm containing conditional branches and loops. As a matter
of fact, there are a lot of functions appearing in scientific computation
problems, which cannot be expressed in a closed mathematical form.

The basic idea behind the \AD of programs is 
very simple. Applying the elementary rules of differentiation to each
statement in the computation-sequence of $f$ yields a new sequence
for the computation of the desired derivatives. The key to this method
is the repeated use of the chain rule \cite{Kede80}.\\[0.3cm]
Why do we use this \AD for computing various 
derivatives of a function $f$\/?
\begin{itemize}
\item the numerical evaluation of the generated computation-sequences is
      much more accurate and faster than the methods based on {\sl Divided
      Differences} (see \cite{Grie89} for a comparison of these methods)
\item the representation of the function and its derivatives as a
      computation-sequence is in general far more compact than a formula
      or even a DAG (Directed Acyclic Graph).
\item we'll be able to handle a larger class of functions by using
      \AD
\item P.~Wolfe \cite{Wolf82} showed that the cost of evaluating a gradient 
      with $n$ components is uniformly bounded and doesn't depend on the number
      $n$ of independent variables.\\
      Hence \AD is really a powerful tool for generating effective code for
      gradient--calculations
\end{itemize}
The disadvantages of a purely symbolical approach are as follows. 
Functions containing conditional branches and loops can't be
handled that easily. Though the function should be representable
by a closed mathematical formula, the more complex this formula gets,
the more time it takes to compute its derivative, the larger its
derivatives become and hence the more computer memory is needed.

We want to show the advantages of \AD by looking at an example.\\
Consider the function $f$ defined by $f_0 = 0 , f_1 = x$
and $f_n = \ln(f_{n-1}+f_{n-2})$ for $n \geq 2$\\
%\begin{small}\begin{verbatim}%that works!!!
%> f5 := ln(ln(ln(x)+ln(x+ln(x)))+ln(ln(x+ln(x))+ln(ln(x)+ln(x+ln(x))))):
%> lprint(diff(f5,x));
%((1/x+(1+1/x)/(x+ln(x)))/(ln(x)+ln(x+ln(x)))+((1
%+1/x)/(x+ln(x))+(1/x+(1+1/x)/(x+ln(x)))/(ln(x)+ln(
%x+ln(x))))/(ln(x+ln(x))+ln(ln(x)+ln(x+ln(x)))))/(
%ln(ln(x)+ln(x+ln(x)))+ln(ln(x+ln(x))+ln(ln(x)+ln(x
%+ln(x)))))
%\end{verbatim}\end{small}
%\begin{small}\begin{verbatim}%that doesn't work!!!
%   > f5 := ln(ln(ln(x)+ln(x+ln(x)))+ln(ln(x+ln(x))+ln(ln(x)+ln(x+ln(x))))):
%   > lprint(diff(f5,x));
%   ((1/x+(1+1/x)/(x+ln(x)))/(ln(x)+ln(x+ln(x)))+((1+1/x)/(x+
%   ln(x))+(1/x+(1+1/x)/(x+ln(x)))/(ln(x)+ln(x+ln(x))))/(ln(x+ln
%   (x))+ln(ln(x)+ln(x+ln(x)))))/(ln(ln(x)+ln(x+ln(x)))+ln(ln(x+
%   ln(x))+ln(ln(x)+ln(x+ln(x)))))
%\end{verbatim}\end{small}
\begin{small}\noindent
\mbox{}\verb!     > f5 := ln(ln(ln(x)+ln(x+ln(x)))+!
\verb!ln(ln(x+ln(x))+ln(ln(x)+!
\verb!ln(x+ln(x))))):!\\
\mbox{}\verb!     > lprint(diff(f5,x));!\\
\mbox{}\verb!     ((1/x+(1+1/x)/(x+ln(x)))/(ln(x)+!
\verb!ln(x+ln(x)))+((1+1/x)/(x+ln(x))+!
\verb!(1/x+(1+1/x)/!\\
\mbox{}\verb!     (x+ln(x)))/(ln(x)+ln(x+ln(x))))/!
\verb!(ln(x+ln(x))+ln(ln(x)+!
\verb!ln(x+ln(x)))))/(ln(ln(x)+!\\
\mbox{}\verb!     ln(x+ln(x)))+ln(ln(x+ln(x))+!
\verb!ln(ln(x)+ln(x+ln(x)))))!\\
\end{small}

This is the representation of \verb/f/ and its derivative by a formula.
Systems that operate on these formulas are very ineffective because the
size of the formula for $f$ and its derivatives grows very large as
a function of $n$.  We need a better representation for $f$.

\noindent Another way to represent $f$ is by a {\bf DAG} ({\bf D}irected 
{\bf A}cyclic {\bf G}raph) which looks like this:
\begin{center} {\large DAG without any sharing}\\
\begin{Large}
\unitlength0.5cm
\begin{picture}(17,15.5) \thicklines
\put(-0.1,0){\makebox(1.2,1.2){\bf x}} 
\put(1.9,0){\makebox(1.2,1.2){\bf x}} 
\put(3.9,0){\makebox(1.2,1.2){\bf x}} 
\put(6.9,0){\makebox(1.2,1.2){\bf x}} 
\put(8.9,0){\makebox(1.2,1.2){\bf x}} 
\put(11.9,0){\makebox(1.2,1.2){\bf x}} 
\put(13.9,0){\makebox(1.2,1.2){\bf x}} 
\put(15.9,0){\makebox(1.2,1.2){\bf x}} 

\put(0.5,1){\vector(0,1){1}}
\put(2.5,1){\vector(0,1){1}}
\put(4.5,1){\vector(0,1){1}}
\put(7.5,1){\vector(0,1){1}}
\put(9.5,1){\vector(0,1){1}}
\put(12.5,1){\vector(0,1){1}}
\put(14.5,1){\vector(0,1){1}}
\put(16.5,1){\vector(0,1){1}}

\put(0,2){\framebox(1,1){ln}}
\put(2,2){\framebox(1,1){ln}}
\put(7,2){\framebox(1,1){ln}}
\put(14,2){\framebox(1,1){ln}}
\put(16,2){\framebox(1,1){ln}}
\put(4,2){\framebox(1,1){+}}
\put(9,2){\framebox(1,1){+}}
\put(12,2){\framebox(1,1){+}}

\put(3,2.5){\vector(1,0){1}}
\put(8,2.5){\vector(1,0){1}}
\put(14,2.5){\vector(-1,0){1}}
\put(4.5,3){\vector(0,1){1}}
\put(9.5,3){\vector(0,1){1}}
\put(12.5,3){\vector(0,1){1}}

\put(4,4){\framebox(1,1){ln}}
\put(9,4){\framebox(1,1){ln}}
\put(12,4){\framebox(1,1){ln}}

\put(2,6){\framebox(1,1){+}}
\put(14,6){\framebox(1,1){+}}

\put(4,5){\vector(-1,1){1}}
\put(1,3){\vector(1,3){1}}

\put(14,7){\vector(-1,1){1}}
\put(3,7){\vector(1,1){1}}

\put(13,5){\vector(1,1){1}}
\put(16,3){\vector(-1,3){1}}

\put(4,8){\framebox(1,1){ln}}
\put(12,8){\framebox(1,1){ln}}

\put(9,10){\framebox(1,1){+}}
\put(9.5,5){\vector(0,1){5}}
\put(12,9){\vector(-2,1){2}}

\put(9,12){\framebox(1,1){ln}}
\put(9.5,11){\vector(0,1){1}}

\put(7,12){\framebox(1,1){+}}
\put(9,12.5){\vector(-1,0){1}}

\put(5,9){\vector(2,3){2}}

\put(7,14){\framebox(1,1){ln}}
\put(8.2,14){\makebox(0.8,1.2){$f_5$}}

\put(7.5,13){\vector(0,1){1}}


\end{picture}
\end{Large}
\end{center}


\noindent
If repeated sub-expressions are shared, the DAG looks like this

\begin{center}

{\large DAG with common subexpressions computed and stored only once}

\begin{Large}
\unitlength0.5cm
\begin{picture}(16,5.5) \thicklines
\put(-1.8,0.2){\makebox(1,1){\bf x}} 
\put(2,0){\framebox(1,1){ln}}
\put(2,2){\framebox(1,1){+}}
\put(2.5,1){\vector(0,1){1}}
\put(3,2.5){\vector(1,0){1}}
\put(3,1){\vector(3,1){3}}

\put(4,2){\framebox(1,1){ln}}
\put(6,2){\framebox(1,1){+}}
\put(5,2.5){\vector(1,0){1}}
\put(5,3){\vector(3,1){3}}
\put(7,2.5){\vector(1,0){1}}

\put(8,2){\framebox(1,1){ln}}
\put(8,4){\framebox(1,1){+}}
\put(8.5,3){\vector(0,1){1}}
\put(9,3){\vector(3,1){3}}
\put(9,4.5){\vector(1,0){1}}

\put(10,4){\framebox(1,1){ln}}
\put(12,4){\framebox(1,1){+}}
\put(11,4.5){\vector(1,0){1}}
\put(13,4.5){\vector(1,0){1}}

\put(14,4){\framebox(1,1){ln}}

\put(-1,1){\vector(3,1){3}}
\put(-1,0.5){\vector(1,0){3}}

\put(15.2,4){\makebox(0.8,1.2){$f_5$}}

\end{picture}
\end{Large}
\end{center}
The Maple system always optimizes the storage of DAG's in this way.
It's obvious that this representation is far more compact.
Hence a differentiation based on this shared DAG will be faster,
exponentially faster in fact.
The running time of the Maple \verb+diff+ function is proportional
to the size of the optimized DAG.

\noindent
But an even more compact representation for $f$, independent of $n$, is to
represent $f$ as a program.  For example, the corresponding Maple--program
for $f$ is
\begin{small}\begin{verbatim}
 > f := proc(x,n) local a,b,h,i;
            a := 0; b := x; for i to n do  h := b; b := ln(a+b); a := h od; b
        end:
\end{verbatim}
\end{small}

\noindent
Computing the derivative of $f$ as a program we obtain

\begin{small}\begin{verbatim}
 > GRADIENT(f,[x]);

  proc(x,n) local a,b,da,db,dh,h,i;
      da := 0; a := 0;
      db := 1; b := x;
      for i to n do
         dh := db; h := b; db := (da+db)/(a+b); b := ln(a+b); da := dh; a := h
      od;
      db
  end
\end{verbatim}
\end{small}

\noindent
The resulting Maple procedure is quite compact even for an arbitrary number $n$!
The cost of evaluating the function $f'$ and the complexity of the 
corresponding nonoptimized DAG will increase by an 
exponential factor if $n$ increases. Whereas the complexity and the
length of the program--representation won't change!\\[0.2cm]
{\large\bf The most effective approach to differentiate 
arbitrary complex functions is the use of Algorithmic Differentiation.}\\[0.2cm]
The first representation of a function by a DAG goes back to Kantarovich
\cite{Kant57}.
Hence this DAG is called a {\sl Kantarovich--Graph}\/! There are two 
possible ways to traverse the graph. One method is obvious. Start at the 
bottom and go through all nodes to the top of the computational--tree.
That's nothing else than passing the corresponding computation--sequence
top--down as usual. Based on this method, the so called 
{\sl forward--mode} or {\sl forward--sweep} of \AD was developed
\cite{Kede80,Grie89,Grie91}.

In contrast to this {\sl forward--mode} the {\sl reverse--mode} or
{\sl reverse--sweep} was first used by Speelpenning \cite{Spee80}. He passed the 
underlying graph top--down and propagated the gradient backwards. He 
also proved that this {\sl reverse--mode} generates automatically piece-wise
optimal computation--sequences (see also \cite{Grie89,Grie91}) in the sence that
all common subexpressions appeared by differentiation are optimized but
not the whole computation--sequence.\\[0.1cm]

In the next section we will describe the {\sl forward--mode} and the implementation
in the Computer Algebra system Maple. Section 3 shows how the function
\verb/GRADIENT/ handles local procedures. In section 4, we want to look at
the syntax of \verb/GRADIENT/. We summarize the class of
Maple--procedures handled by it and the three subsections describe the 
different functions based on \verb/GRADIENT/. In section 5 we draw
the conclusion from this approach in Maple.

% ////////////////////////////////////////////////////////////////////////

\section{The forward-mode in Maple}

The function \verb/GRADIENT/ presented in this paper is based on a very first
version of the function \verb/PD/ \cite{Mona91}.
Given a function $f$ depending on $\vec{x}=(x_1,\ldots,x_n)$ represented
by this program\\
\hspace*{1.2cm}
\parbox{7cm}
{
    \begin{tabbing}
    {\bf Original Program}\\[0.2cm]
    {\bf begin}\\
    \qquad \= {\sl statement}$_{\,1}${\bf ;}\\
    \> {\sl statement}$_{\,2}${\bf ;}\\
    \> \qquad $\vdots$\\
    \> {\sl statement}$_{\,m}$\\
    {\bf end:}
    \end{tabbing}
}\hfill
\parbox{1cm}{\begin{eqnarray} \label{MyOrigProg} \end{eqnarray}}\\
With the following notations:
\begin{itemize}
\item $\vec{x}=(x_1,\ldots,x_n)$ are the independent variables
\item $T$ is a set of all necessary temporary variables which depend on 
          $\vec{x}$ 
\item $B$ is a set of the union of all mathematical--functions, the elementary
          arithmetic operations $+,-,\cdot,/$ and all possible
          compositions of them
\item $C$ is the set of all constants
\item {\sl statement}$_i\/$ is a single program--statement which corresponds
          to one of the five types listed below:\\
\parbox{8cm}{
   \begin{itemize}
   \label{TypList}
   \item[i)] $t := h(C,T,\vec{x})$; \qquad where \quad $t \in T,h \in B$
   \item[ii)] \verb/if/--control--statement
   \item[iii)] \verb/for/--statement
   \item[iv)] \verb/RETURN/--statement 
   \item[v)] a single mathematical expression
   \end{itemize}}\hfill
\parbox{1cm}{\begin{eqnarray} \label{condition} \end{eqnarray}}\\
\end{itemize}
\parbox[b]{1cm}{\bf Note:}
\parbox[t]{14cm}{The structure of the different branches of an 
                \verb/if/--statement and the body of a \verb/for/--loop 
                is also a statement--sequence like (\ref{MyOrigProg}).}\\[0.2cm]
The source--code transformation of the original program in order to
get a program for the computation of the gradient of $f$ using the
{\sl forward--mode} isn't very difficult.

The sign $\nabla$ denotes the gradient with respect to all variables
$x_1,\ldots,x_n$. While we're passing the Kantarovich--Graph of $f$ 
bottom--up, we are adding new nodes to it such that we get a larger
graph representing the calculation of the gradient of $f$.
From this graph we can derive the\\
\hspace*{1.2cm}
\parbox{7cm}
{
  {\bf Differentiated Program}\\
  \hspace*{0.22cm}
   \parbox{5cm}{
    \begin{tabbing}
    {\bf begin}\\
    \qquad \= {\bf for} \= $i$ {\bf from} $1$ {\bf to} $n$ {\bf do}\\
     \> \> $\nabla x_i = e_i$\\
      \> {\bf od;}\\
       \> $\nabla\!${\sl statement}$_{\,1}${\bf ;}\\
        \> $\nabla\!${\sl statement}$_{\,2}${\bf ;}\\
         \> \qquad $\vdots$\\
          \> $\nabla\!${\sl statement}$_{\,m}$\\
           {\bf end:}
           \end{tabbing}}
}\hfill
\parbox{1cm}{\begin{eqnarray}\label{MyDiffProg}\end{eqnarray}}\\
Where $e_i$ denotes the $i$-th Cartesian basis vector in $\rR^n$. 
The set $G$ shall contain the gradients of all variables of $T$. 
All $\nabla\!${\sl statement}$_{\,1}$,\ldots,$\nabla\!${\sl 
statement}$_{\,m}$ shall be the ``gradients'' of the corresponding
statements. Each {\sl statement}$_{\,i}\:$ will be replaced by $\nabla\!${\sl 
statement}$_{\,i}\:$ which is one of the following five types:
\begin{itemize}
\item[i)] $\nabla t :=\sum\limits_{z\in\{x_1,\ldots,x_n\}\cup T}
                       \frac{\partial h(C,T,\vec{x})}
                      {\partial z} \cdot\nabla z$;\\
                 $t := h(C,T,\vec{x})$; \qquad 
                 \qquad \qquad \qquad where $t\in T,\nabla t\in G,h\in B$
\item[ii)] \verb/if/--control--statement where the statement--sequences
           in all branches were replaced by their "gradients"
\item[iii)] \verb/for/--loop where the iterated statement--sequence was 
            replaced by its ``gradient''
\item[iv)] \verb/RETURN/--statement which returns the calculated gradient
\item[v)] \verb/RETURN/--statement which returns the gradient of the 
          corresponding mathematical expression
\end{itemize}
\parbox[b]{1cm}{\bf Note:}
\parbox[t]{14cm}{We have to compute $\nabla t$ {\bf before} the computation of $t$ if
in the assignment $t := h(C,T,\vec{x})$, $h(C,T,\vec{x})$ itself 
depends on $t$.}\\[0.3cm]

The first thing we had to do was to choose a suitable data structure
to represent a gradient. We decided to use the array data--type because it's 
very close to the mathematical formalism and allows you to build a Jacobian-- or
even a Hessian--procedure basing on \verb/GRADIENT/ quite easily. 
Another advantage is that the user can include the generated 
``gradient--procedures'' in his own programs because  their return types
are always well defined.\\[0.1cm]

Assume that $f=f(x_1,\ldots,x_n) : \rR^n \to \rR^m$. The corresponding procedure
\verb/f/ shall return a one-dimensional array with $m$ entries. By using
\verb/GRADIENT/ we get a new procedure for $\nabla\! f$ that returns
the Jacobian matrix of $f$ with $m \cdot n$ entries.

Lets look at two simple examples to illustrate the use of arrays:\\[0.2cm]
{\bf Example 1}:
\begin{small}\begin{verbatim}
> f := proc(x,y) local a; a := array(1 .. 2);
           a[1] := x^2*y^2;  a[2] := x^2+y^2;  a
       end:
> eval(f(x,y));
                                  2  2   2    2
                               [ x  y , x  + y  ]
> df := GRADIENT(f);

df := proc(x,y) local da;
          da := array(1 .. 2,1 .. 2);
          da[1,1] := 2*x*y^2;
          da[1,2] := 2*x^2*y;
          da[2,1] := 2*x;
          da[2,2] := 2*y;
          da
      end
> eval(df(x,y));
                               [      2     2   ]
                               [ 2 x y   2 x  y ]
                               [                ]
                               [   2 x     2 y  ]
\end{verbatim}\end{small}
{\bf Example 2:}
\begin{small}\begin{verbatim}
> f := proc(x,y,n) local a,i; a := array(1 .. n); 
           for i to n do a[i] := (x*y)^i od;  a
       end:
> eval(f(x,y,2));
                                         2  2
                                 [ x y, x  y  ]
> df := GRADIENT(f,[x,y]);

df := proc(x,y,n) local a,i,da;
          da := array(1 .. n,1 .. 2);
          for i to n do
              da[i,1] := (x*y)^i*i/x;  # diff(a[i],x)
              da[i,2] := (x*y)^i*i/y;  # diff(a[i],y)
              a[i] := (x*y)^i
          od;
          da
      end

> eval(df(x,y,2));
                               [    y       x   ]
                               [                ]
                               [      2     2   ]
                               [ 2 x y   2 x  y ]
\end{verbatim}\end{small}
The second example shows that we even don't have to fix
the arraybounds! If we verify the printed results we'll see that these are just
the Jacobian matrices of the underlying functions! By an extensive use of this
array data--structure we don't have to care about multivalued functions. We
just call the \verb/GRADIENT/--function and it will create automatically the
correct Jacobian matrix of \verb/f/.\\[0.3cm]
If we want to compute the gradient of an arbitrary mathematical expression the
only thing we have to do is to differentiate the function $h(C,T,\vec{x})$
where $h \in B$ with respect to the desired variables.\\[0.1cm]
The result is a vector of the form $\sum\limits_{z\in\{x_1,\ldots,x_n\}\cup T}
  \frac{\partial h}{\partial z}\cdot \nabla z$.\\[0.1cm]
Lets look at the following example to explain the procedure.\\[0.2cm]
{\bf Example}:
\[ h=h(x,y,a,b):=\sin(x)y+a^{2}b \]
\begin{quote}
-- $x,y$ are the independent variables\\
-- $a,b$ are temporary variables depending on $x$ and $y$\\[0.1cm]
Assume that $da$ and $db$ are declared as one--dimensional arrays with two
entries denoting the partial derivatives of $a$ and $b$.\\
If we calculate the {\it Total Derivative} of $h$ we get: 
\[ \nabla h = \frac{\partial h}{\partial x}\nabla x +
              \frac{\partial h}{\partial y}\nabla y +
              \frac{\partial h}{\partial a}\nabla a +
              \frac{\partial h}{\partial b}\nabla b 
\]
\begin{eqnarray*}
 \nabla h & = & \cos(x)y\cdot\left(\! \begin{array}{c} 1 \\ 0
\end{array}\!\right)+\sin(x)\cdot\left(\!
\begin{array}{c} 0 \\ 1 \end{array}\!\right) + 2\,a\,b\cdot\left(\!
\begin{array}{c} \frac{\partial a}{\partial x} \\ 
 \frac{\partial a}{\partial y} \end{array}\!\right) + a^{2}\cdot\left(\!
\begin{array}{c} \frac{\partial b}{\partial x} \\ \frac{\partial b}{\partial y}
 \end{array}\!\right) 
\\
& = & \left(\begin{array}{c} \cos(x)y+2\,a\,b\cdot
              \mbox{{\sl da}[1]}+a^{2}\cdot\mbox{{\sl db}[1]} \\
              \sin(x)+2\,a\,b\cdot\mbox{{\sl da}[2]}
               +a^{2}\cdot\mbox{{\sl db}[2]}
            \end{array}
      \right)
\end{eqnarray*}
respectively:
\end{quote}
The following example shall illustrate the differentiation of
various kinds of statements supported by \verb/GRADIENT/:
\begin{small}\begin{verbatim}
 > f := proc(x,y,n) local i,t;
            t := x*y; for i to n do  t := t*x*y od;
            if x < y then RETURN(t*x,t*y) else RETURN(t,x*y) fi
        end:
 > GRADIENT(f,[x,y]);

  proc(x,y,n) local dt,i,t;
      dt := array(1 .. 2); dt[1] := y; dt[2] := x; t := x*y;
      for i to n do
          dt[1] := dt[1]*x*y+t*y; dt[2] := dt[2]*x*y+t*x; t := t*x*y
      od;
      if x < y then RETURN(dt[1]*x+t,dt[2]*x,dt[1]*y,dt[2]*y+t)
      else RETURN(dt[1],dt[2],y,x)
      fi
  end
\end{verbatim}\end{small}

% ////////////////////////////////////////////////////////////////////////////

\section{Differentiation of local procedures}

The differentiation of local procedures could be done in two different ways.
One is, for reasons of consistency, quite obvious. We call recursively the
function \verb/GRADIENT/ in order to differentiate the local procedure with
respect to only one parameter. This yields a vector (gradient) of
differentiated procedures. But this approach is definitely very time-- and
memory-consuming. A further disadvantage is that for each function call of such
a component--procedure a full evaluation is unavoidable!

The other variant of differentiating local defined functions is that we
immediately compute its gradient, such that we have to do only one function
call to get all the derivatives we want. The run--time of the
\verb/GRADIENT/--function and the generated program will be dramatically
shortened.\\[0.1cm]
\parbox[b]{2.1cm}{\bf Important:}
\parbox[t]{14cm}{Local procedures will be differentiated always with respect to
all their parameters. Hence it's not possible to pass an additional parameter
denoting for example the dimension of an array yet. By using global defined
variables this disadvantage can easily be removed. It is planed to extend
the algorithm such that it should be able to detect the independent parameters
of local procedures. In that case it will be possible to pass further
parameters. Hence the local procedures won't be differentiated with respect
to all their arguments.}\\[0.1cm]

In order to get efficient program sequences the user should take care in
avoiding several calls to the same function with the same arguments. 
For example 
\begin{small}\begin{verbatim}
 > f := proc(x,y) local a,d;
            d := proc(u,v) local b; 
                     b := array(1..2); 
                     b[1] := u^3*v; b[2] := u*v^2; 
                     b 
                 end;
            a := d(x,y);
            a[1]*a[2]*x^2*y
        end:
\end{verbatim}\end{small}

For the function \verb/GRADIENT/ it's important to know the type (data
structure) of the value calculated by a local function. If it doesn't know
this, it'll be impossible to compute the correct total derivative. Consider the
following example
\begin{small}\begin{verbatim}
 > f := proc(x,y) local s,d;
            d := proc(u,v,w) local a; a := u^2*v*w^3; a*u*v end;
            s := d(x,x*y,y);
            s*x^2*y
        end:
 > df := GRADIENT(f);

  df := proc(x,y) local d,dd,ds,lf1,s;
            ds := array(1 .. 2);
            dd := proc(u,v,w) local a,da,grd;
                      grd := array(1 .. 3);
                      da := array(1 .. 3);
                      da[1] := 2*u*v*w^3;
                      da[2] := u^2*w^3;
                      da[3] := 3*u^2*v*w^2;
                      a := u^2*v*w^3;
                      grd[1] := da[1]*u*v+a*v;
                      grd[2] := da[2]*u*v+a*u;
                      grd[3] := da[3]*u*v;
                      grd
                  end                         ;
            d := proc(u,v,w) local a; a := u^2*v*w^3; a*u*v end;
            lf1 := dd(x,x*y,y);
            ds[1] := lf1[1]+lf1[2]*y;
            ds[2] := lf1[2]*x+lf1[3];
            s := d(x,x*y,y);
            RETURN(ds[1]*x^2*y+2*s*x*y,ds[2]*x^2*y+s*x^2)
        end
\end{verbatim}\end{small}
%
%> f(x,y);df(x,y);
%
%                                      7  6
%                                     x  y
%
%                                   6  6     7  5
%                                7 x  y , 6 x  y

The local function \verb/d/ has three parameters and returns a single value.
The differentiated function \verb/dd/ returns now the gradient as an array with
three entries. This gradient is stored in the temporary variable \verb/lf1/. In
the original program the variable \verb/s/ holds the value \verb/g(x,y,y)/,
\verb/ds/ should denote now the gradient of \verb/g(x,y,y)/ with respect to
\verb/x/ and \verb/y/. The componentwise assignments to \verb/ds[1]/ and 
\verb/ds[2]/ lead to the correct result.

% ///////////////////////////////////////////////////////////////////////////////

\section{The GRADIENT procedure}

The Maple--procedure \verb/GRADIENT/ computes $\nabla f$ where $f$ is input as a
Maple--procedure.
It is possible to call the \verb/GRADIENT/ procedure in three different ways.
\begin{itemize}
\item[i)] \verb/GRADIENT(f)/ 
\item[ii)] \verb/GRADIENT(f,p)/ 
\item[iii)] \verb/GRADIENT(f,[x1,..,xn])/ 
\end{itemize}
With the following definitions:\\[0.1cm]
\parbox[b]{2cm}{\quad \verb/f/ :}
\parbox[t]{14cm}{a Maple procedure which computes $f$}\\
\parbox[b]{2cm}{\quad \verb/p/ :}
\parbox[t]{14cm}{a positive integer denoting that the procedure 
                  \verb/f/ will be differentiated with respect to the first
                  \verb/p/ parameters in the parameterlist of
                  \verb/f/.  Other parameters will be interpreted as
                  constants.}\\[0.1cm]
\parbox[b]{2cm}{\verb/ x1,..,xn/ :}
\parbox[t]{14cm}{a list of variables with respect to which the gradient
                         will be computed}\\[0.2cm]
{\bf \underline{Restrictions to the Maple procedures}}
\begin{itemize}
\item procedure--options as \verb/remember/, \verb/trace/ etc. aren't allowed
\item it's not allowed to use \verb/args/, \verb/nargs/ nor quotes
\item recursive procedures are not supported yet
\item assignments to global variables are forbidden too
\item loop variables have to be constants
\item only the \verb/array/ data structure is supported
\item array definitions mussn't contain options, index functions nor
      initialization lists
\end{itemize}
{\bf \underline{What is allowed in Maple procedures}}
\begin{itemize}
\item operations with the \verb/array/ data type
\item assignment statements
\item \verb/for/--loops, \verb/while/--loops and \verb/if/--statements
\item the statements \verb/RETURN/, \verb/ERROR/, \verb/break/, \verb/print/ and
      \verb/lprint/ 
\item the definition of local procedures
\item the use of the elementary arithmetic operations \verb/+, -, *,/
      \verb+/, ^+ 
\item any mathematical function Maple knows how to differentiate
\end{itemize}
\vspace{0.3cm}

\subsection{The JACOBIAN procedure}

The Jacobian matrix $J$ of the function 
\begin{equation} \label{VektorFunc}
 \begin{array}{cccc} 
  \begin{array}{c} \vec{f} : \\ \ \end{array}
   \begin{array}{c} \rR^n \\ \vec{x} \end{array}
    \begin{array}{c} \to \\ \mapsto \end{array}
     \begin{array}{c} \rR^m \\ \vec{y}=\vec{f}(\vec{x}) \end{array}
      \end{array}
\end{equation}
is defined by
\begin{equation} \label{JacTerme}
  J_{i,j} = \frac{\partial y_i}{\partial x_j} \qquad \mbox{where}
  \qquad \begin{array}{c} 1\leq i \leq m \\ 1\leq j \leq n
          \end{array}
   \end{equation}
These are naturally the entries in the Jacobian matrix $J_{\vec{f}}$ of the
function $\vec{f}$.\\[0.2cm]
Given a program that computes the function $\vec{f}$. Generate with the aid of
\AD a new program for the simultaneous computation of the terms
(\ref{JacTerme}).

We've already shown that it's possible to generate automatically the desired
Jacobian matrix by using arrays.
In many scientific applications it's often the case that each component of
$\vec{f}$ is given by an own function $f_i$ $(1 \leq i \leq m)$. Thus we have
$m$ functions $f_1,\ldots,f_m$ depending all on $\vec{x}=(x_1,\ldots,x_n)$ and
each represented by a program.
The only thing we have to do is to differentiate each function $f_i$ and
collect the results as local functions in a new procedure.

The implemented function \verb/JACOBIAN/ is based on \verb/GRADIENT/. Thus each of
the procedures \verb/f1,..,fm/ must satisfy the conditions (\ref{condition}).
The syntax of the \verb/JACOBIAN/ procedure is similar to the \verb/GRADIENT/
procedure
\begin{itemize}
\item[i)]   \verb/JACOBIAN([f1,..,fm])/
\item[ii)]  \verb/JACOBIAN([f1,..,fm],p)/
\item[iii)] \verb/JACOBIAN([f1,..,fm],[x1,..xn])/
\end{itemize}
Where \verb/f1,..,fm/ are Maple--procedures.
\subsubsection*{Example:}
\begin{small}\begin{verbatim}
 > f1 := proc(x,y) local a; a := x^3*y^2; a*x*y end:
 > f2 := proc(u,v) local a; a := 3*u^2*v^3; 4*a*v^2 end:
 > jacf := JACOBI([f1,f2]);

 jacf := proc(t1,t2) local gf,jm,i,h;
            jm := array(1 .. 2,1 .. 2);
            gf := array(1 .. 2);
            gf[1] := proc(x,y)
                     local a,da,grd;
                         grd := array(1 .. 2);
                         da[1] := 3*x^2*y^2;
                         da[2] := 2*x^3*y;
                         a := x^3*y^2;
                         grd[1] := da[1]*x*y+a*y;
                         grd[2] := da[2]*x*y+a*x;
                         grd
                     end                         ;
            h := gf[1](t1,t2);
            for i to 2 do  jm[1,i] := h[i] od;
            gf[2] := proc(u,v)
                     local a,da,grd;
                         grd := array(1 .. 2);
                         da[1] := 6*u*v^3;
                         da[2] := 9*u^2*v^2;
                         a := 3*u^2*v^3;
                         grd[1] := 4*da[1]*v^2;
                         grd[2] := 4*da[2]*v^2+8*a*v;
                         grd
                     end                             ;
            h := gf[2](t1,t2);
            for i to 2 do  jm[2,i] := h[i] od;
            jm
        end

 > f1(a,b),f2(a,b);
                                 4  3      2  5
                                a  b , 12 a  b
 > eval(jacf(a,b));

                             [    3  3      4  2 ]
                             [ 4 a  b    3 a  b  ]
                             [                   ]
                             [       5      2  4 ]
                             [ 24 a b   60 a  b  ]
\end{verbatim}\end{small}

\subsection{The HESSIAN procedure}

The Hessian matrix $H$ of a function $f : \rR^n \to \rR$ is the matrix
\begin{eqnarray} \label{Hessterms}
  H_{i,j} = \frac{\partial^2 f}{\partial x_i\partial x_j}
\quad\mbox{ where }\quad  1 \leq i,j \leq n
\end{eqnarray}
Given a program that computes the function $f$. Generate with the aid of \AD a
new program for the simultaneous computation of the terms
(\ref{Hessterms}).

The implemented function \verb/HESSIAN/ is based on \verb/GRADIENT/ too. For that
reason the procedure \verb/f/ has to satisfy the conditions (\ref{condition}).
The syntax of the \verb/HESSIAN/ procedure is exactly the same as it is for the
\verb/GRADIENT/ procedure.
\begin{itemize}
\item[i)] \verb/HESSIAN(f)/ 
\item[ii)] \verb/HESSIAN(f,p)/ 
\item[iii)] \verb/HESSIAN(f,[x1,..,xn])/ 
\end{itemize}

\subsubsection*{Example:}
\begin{small} \begin{verbatim}
> f := proc(x,y) local a; a := x^2*y; a^3*y end:
> f(x,y);
                                     6  4
                                    x  y

>  hessf := HESSIAN(f);
 
 hessf := proc(x,y) local a,t1,da,dda;
              da := array(1 .. 2);
              dda := array(1 .. 2,1 .. 2);
              dda[1,1] := 2*y;
              dda[1,2] := 2*x;
              da[1] := 2*x*y;
              da[2] := x^2;
              a := x^2*y;
              t1 := array(1 .. 2,1 .. 2);
              t1[1,1] := 6*a*y*da[1]^2+3*a^2*y*dda[1,1];
              t1[1,2] := 6*a*y*da[2]*da[1]+3*a^2*y*dda[1,2]
                         +3*a^2*da[1];
              t1[2,1] := t1[1,2];
              t1[2,2] := 6*a*y*da[2]^2+6*a^2*da[2];
              t1
          end

 > eval(hessf(x,y));
                             [     4  4      5  3 ]
                             [ 30 x  y   24 x  y  ]
                             [                    ]
                             [     5  3      6  2 ]
                             [ 24 x  y   12 x  y  ]
\end{verbatim}\end{small}

\subsection{The TAYLOR procedure}

Lots of scientific problems are solved by approximations using Taylor series.
The function \verb/TAYLOR/ generates a procedure for the computation of all
terms of the following form
\begin{eqnarray} \label{Taylorterms}
\frac{\partial^k f}{\partial^{k_1}x_1\cdot\partial^{k_2}x_2
\cdots\partial^{k_n}x_n}
\quad\mbox{ where }\quad
\begin{array}{c} k_i \in {I\!\!N}_0 \quad\mbox{with}\quad i=1,\ldots,n \\
                 k := \sum\limits_{i=1}^n k_i \leq d
\end{array}
\end{eqnarray}
Note that \verb/TAYLOR/ takes care in the order of differentiation and
eliminates all redundant derivatives!\\[0.2cm]
The procedure \verb/f/ representing the function $f$ has to satisfy as well the
conditions (\ref{condition}).

The syntax of the \verb/TAYLOR/ procedure is similar to the \verb/GRADIENT/
procedure
\begin{itemize}
\item[i)]   \verb/TAYLOR(f,d)/
\item[ii)]  \verb/TAYLOR(f,p,d)/
\item[iii)] \verb/TAYLOR(f,[x1,..,xn],d)/
\end{itemize}
The parameter \verb/d/ denotes the maximal total degree of the derivatives in
(\ref{Taylorterms}).

\subsubsection*{Example:}
\begin{small}\begin{verbatim}
 > f := proc(x,y) local a; a := x^2*y; a^3*y end:
 > tayf := TAYLOR(f,2);
 tayf := proc(x,y) local a,t1,t2,t3,da,dda;
     da := array(1 .. 2);
     dda := array(1 .. 2,1 .. 2);
     dda[1,1] := 2*y;
     dda[1,2] := 2*x;
     da[1] := 2*x*y;
     da[2] := x^2;
     a := x^2*y;
     t1 := array(1 .. 3);
     t1[1] := a^3*y;
     t2 := array(1 .. 2);
     t2[1] := 3*a^2*y*da[1];
     t2[2] := 3*a^2*y*da[2]+a^3;
     t1[2] := t2;
     t3 := array(1 .. 3);
     t3[1] := 6*a*y*da[1]^2+3*a^2*y*dda[1,1];
     t3[2] := 6*a*y*da[2]*da[1]+3*a^2*y*dda[1,2]
              +3*a^2*da[1];
     t3[3] := 6*a*y*da[2]^2+6*a^2*da[2];
     t1[3] := t3;
     t1
 end

 > eval(tayf(x,y));
                               6  4
                            [ x  y , t2, t3 ]

 > eval("[2]);
                               5  4     6  3
                          [ 6 x  y , 4 x  y  ]

 > eval(""[3]);
                           4  4      5  3      6  2
                     [ 30 x  y , 24 x  y , 12 x  y  ]

\end{verbatim}\end{small}

\section{Conclusion}

As you can see the use of a Computer Algebra system for doing \AD is a
good alternative to an approach in a object oriented programming environment
(see \cite{Rall84,ADOL91}). You don't have to worry about differentiation and 
simplification of mathematical expressions. That's all done by the system.
You can put all of your effort in extending your algorithm such that it handles
various data types and high order derivatives. It is possible to produce Fortran-
or even C- code of Maple--procedures, such that the user may program his function
in Maple, then uses the AD--algorithm to differentiate it and afterwards convert
the generated procedure to C in order to implement it in other C--routines. The
user can be assured that his routines for the computation of $f$ and even the one
for the computation of its derivatives are optimized and very efficient codes.
Another advantage is the one that the procedures generated by \AD can
immediately be used for further computations in Maple.

The most difficult part is the optimization of the resulting computation--sequence
by using the ``forward--mode'' of \AD. The implemented \verb/GRADIENT/ 
routine does constant-folding and a few optimizations to eliminate unnecessary statements.
Some mathematical simplifications are done by Maple automatically.
But this is not sufficient. In the near future this optimizations should be 
done by Maple's \verb/optimize/ routine.

The next step to be done is the implementation of the ``reverse--mode'' of \AD.
The implementation of this ``reverse--mode'' is much harder than the one for the 
``forward--mode'' because the whole computational tree of $f$ must be stored
such that it can be passed backwards. That means that we have to take special 
care in \verb/if/--statements and \verb/for/--loops.

Furthermore it would be useful if \verb/GRADIENT/ handled recursive procedures.
This feature has not been implemented yet because it requires an additional
pass through the original procedure.

\begin{thebibliography}{99}
\label{Lit}
  \bibitem{Beda59} L.~M.~Beda et al, {\sl Programs for Automatic 
                   Differentiation for the Macine BESM}, Inst. Precise 
                   Mechanics and Computation 
                   Techniques, Academy of Science, Moscow, 1959
  \bibitem{Mapl91} B.~Char, K.~Geddes, G.~Gonnet, B.~Leong, M.~Monagan, 
                   S.~Watt, {\sl Maple V Language Reference Manual}, 
                   Springer, 1991
  \bibitem{Grie89} A.~Griewank, {\sl On Automatic Differentiation}, 
                   Mathematical Programming, edited by M.~Iri and K.~Tanabe, 
                   Kluwer Academic Publishers, pp. 83-107, 1989
  \bibitem{Grie91} A.~Griewank, {\sl The Cain Rule Revisited in Scientific 
                   Computing}, SIAM News May 1991, part I pp. 20-21, SIAM 
                   News July 1991, part II pp. 8-9, p. 24
  \bibitem{ADOL91} A.~Griewank, D.~Juedes, J.~Srinivasan, C.~Tyner,
                   {\sl ADOL--C, A package for the automatic differentiation
                   of algorithms written in C/C++}, ACM Trans.~Math.~Software,
                   to appear
  \bibitem{Iri84}  M.~Iri, {\sl Simultaneous Computation of Functions, 
                   Partial Derivatives and Estimates of Rounding Errors -- 
                   Complexity and Practicality}, Japan Journal of Applied 
                   Mathematics, Vol.~1, No.~2, pp. 223-252, 1984
  \bibitem{Joss76} J.~Joss, {\sl Algorithmisches Differenzieren}, 
                   Diss. ETH 5757, 1976
  \bibitem{Jued91} D.~W.~Juedes, {\sl A Taxonomy of Automatic 
                   Differentiation Tools}, in A.~Griewank and G.F.~Corliss, 
                   editors, {\sl Automatic Differentiation
                   of algorithms: Theory, Implementation and Application},
                   SIAM, Philadelphia, 1991, 315--329
  \bibitem{Kant57} L.~V.~Kantarovich, {\sl On a mathematical symbolism 
                   convenient for performing machine calculations}, 
                   Dokl. Akad. Nauk SSSR, 
                   Vol.~113, pp. 738-741, 1957
  \bibitem{Kede80} G.~Kedem, {\sl Automatic Differentiation of 
                   Computer Programs},
                   ACM TOMS, Vol.~6, No.~2, pp. 150-165, June 1980
  \bibitem{Mona91} M.~Monagan, {\sl Program PD}, Maple share library, 1991
  \bibitem{Mose74} J.~Moses, {\sl The Evolution of Algebraic Manipulation
                   Algorithms}, IFIP Congress 1974
  \bibitem{Rall81} L.~B.~Rall, {\sl Automatic Differentiation -- 
                   Techniques and Applications}, Springer Lecture Notes 
                   in Computer Science, Vol.~120,1981
  \bibitem{Rall84} L.~B.~Rall, {\sl Differentiation in PASCAL-SC : 
                   Type GRADIENT}, ACM TOMS, Vol.~10, No.~2, pp. 161-184, 
                   June 1984
  \bibitem{SIAM92} A.~Griewank, G.F.~Corliss, {\sl Automatic Differentiation
                   of algorithms: Theory, Implementation and Application},
                   SIAM, Philadelphia, 1991
  \bibitem{Spee80} B.~Speelpenning, {\sl Compiling fast Partial Derivatives of 
                   Functions given by Algorithms}, Ph.D. Dissertation, Department 
                   of Computer Science, University of Illinois at 
                   Urbana-Champaign, Urbana, IL 61801
  \bibitem{Weng64} R.~E.~Wengert, {\sl A simple Automatic Derivative 
                   Evaluation Program}, Com. ACM, Vol.~7, pp. 463-464, 1964
  \bibitem{Wolf82} P.~Wolfe, {\sl Checking the Calculation of 
                   Gradients}, ACM TOMS, Vol.~6, No.~4, pp. 337-343, 1982
\end{thebibliography}

\end{document}

