A wide range of numerical problems can be efficiently solved using derivatives in one way or the other. While from a strictly mathematical point of view the derivative is a well-defined object, its computation is anything but trivial.
A classical approach is finite differences. Let f be differentiable at some point x. Clearly, for a certain h small enough,
The problem with finite differences is twofold. One issue – probably the more important one – is accuracy. Being necessarily an approximation, its quality largely depends on a sensible choice of h. Large values, obviously, make for a poor estimate of the actual derivative; small ones, on the other hand, are prone to computational artefacts such as cancellation. While there are strategies to cope with this dilemma they normally do so – and that is the second concern – at the expense of additional evaluations of your function. Central differences, for instance, requires a total of 2n evaluations, where n is the dimension of the domain space. If, to make matters worse, the computation is carried out within an iterative loop, you forfeit a good deal of the algorithmic efficiency that may have motivated the use of derivatives in the first place.
The concept of Automatic Differentiation is altogether different from the above. Unlike finite differences, it provides a means to analytically compute the derivative of a function at a given inner point of its domain. A straightforward approach – the one implemented by the extension – is to introduce a new data-type, often referred to in literature as differential number or gradient. Basically, this is a compound of the value itself and the associated derivative. The fundamental idea is to define the common operators on the set of differential numbers according to the well-know rules of elementary calculus. Hence, multiplication becomes
and so on for the remaining cases. Now consider that a function, in practice, is implemented by a computer program, which in turn is made up of discrete instructions. Control flow may bifurcate depending on switch-statements, but, no matter how complex its structure, eventually it is a sequence of elementary operations. By overloading all or most of these in the above described manner you create an ideally complete algebra of differential numbers, where
Thus, all you have to do is create an initial gradient and pass it on to the computer program, which will then construct the derivative Df(x) along with the output f(x) simultaneously.
With AD you elude the two principal drawbacks of numerical differentiation outlined previously. First of all, it is more reliable in that you no longer have to worry about approximation errors. Although accuracy, of course, is ultimately bounded by machine precision, it can make a difference if you get 16 instead of, say, 10 correct figures. The other advantage is maybe less apparent and of minor relevance to most users. However, in cases where cost is a non-negligible factor, it may be not indifferent that the number of evaluations does not scale with the problem size. Whether your function depends on 5 or, say, 500 variables, one pass will do either way. Due to the computational overhead implied by every single operation this comes at the price of a slowed-down execution during that single pass. We shall rely on vectorized code for a good performance here.
Today Automatic Differentiation is a widely used technology in both industry and academic sience. Implementations cover almost every language or application commonly used for numerical computations, the most popular being Fortran, C, and Matlab®. For further discussion of the topic and relevant links see, for instance, the INRIA site dedicated to AD.
Copyright© 2006, 2007 Thomas Kasper
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License
for more details.
The maintainer will do his best to retain upward compatibility, though at the current pace of releases this seems an audacious promise.
Latest release: ad-1.0.2.tar.gz (gzipped archive, 175 Kb), a compiled version for windows is available here
Install by typing pkg install ad-x.x.x.tar.gz at the octave prompt. For
advanced options and general information about the new package manager invoke the
online documentation with help pkg or consult the Octave-Forge
website.
It is recommended that you run the integrated testscript after installation to make sure the entire AD-functionality is available to you. Be prepared that it takes a couple of seconds before the results are reported. (Example below for a flawed 2.9.16)
octave:1> fid = fopen ("ad.log", "wt");
octave:2> __ga__ (fid)
PASSES 289 out of 289 tests (7 expected failures)
Note also that the vast majority of tests are statistical and involve randomly generated data. This may occasionally result in noise exceeding the specified tolerance. Do let me know if any of the tests repeatedly fail on your system. Before you report a bug, however, please check the log does not already classify it as a known issue.
D provides an intuitive interface to
AD-functionality while hiding the ugly and potentially confusing details from
the user. Let us have a look at its signature:
octave:3> help D
-- Function File: [Y, J] = D (F, X, VARARGIN)
Evaluate F for a given input X and compute the jacobian, such that
d
J(i,j) = ----- Y(i) where Y = F (X, VARARGIN{:})
dX(j)
If X is complex, the above holds for the directional derivatives
along the real axis
Derivatives are computed analytically via Automatic Differentiation
See also: use_sparse_jacobians.
Note that F must be a function handle, not a character string. A simple
use-case scenario could look as follows. Suppose you want to find a root of the non-linear
function
octave:4> function y = foo (x)
> y(1) = 100 * (x(1) - x(2)^2)^2;
> y(2) = (1 - x(1))^2;
> y = y(:) / 2;
> endfunction
One way to go about this is to iteratively refine a random initial guess by Newton-steps
octave:5> x = rand (2, 1);
octave:6> for k = 1:30, [y, J] = D (@foo, x); x = x - J \ y; endfor
octave:7> x, res = norm (foo (x))
x =
1.0000
1.0000
res = 1.1696e-009
Note that D is a mere convenience function which
wraps up the steps outlined in the introductory section. Thus, [y, J] = D (@F, x)
essentially is a shortcut for
result = F (gradinit (x));
y = result.x;
J = result.J;
gradinit, you specify the independent variables to differentiate with
respect to along with their initial values. In the resulting gradient you find the argument
x augmented by the jacobian which evaluates to the identity matrix of size
numel(x)
octave:5> g = gradinit ([-1; 2])
g =
value =
-1
2
(partial) derivative(s) =
1 0
0 1
Gradients represent a class of their own and are listed as such by the
interpreter. Use isgradient for type-checking:
octave:8> who -long g
*** local user variables:
Prot Name Size Bytes Class
==== ==== ==== ===== =====
rwd g 2x1 48 gradient
Total is 2 elements using 48 bytes
octave:9> isgradient (g)
ans = 1
Each of the two members (value and partial derivatives) can be accessed by
suffixing the variable with .x and .J respectively.
(For obvious reasons, however, they should only be read out and never be assigned
to directly.) Analytical expressions in one or more variables of type gradient
automatically evaluate to gradients:
octave:10> foo (g)
ans =
value =
1250
2
(partial) derivative(s) =
-500 2000
-2 -0
At any time their members satisfy g.J(i,j)
= d/dx(j)[g.x(i)], with x(j) the variables previously passed
to gradinit. Beware that this relation is independent of shape and
extends to arrays of arbitrary dimension. Thus, operations which preserve the
linear order of elements (like reshape, for instance, or a transposal on
column-vectors) do not alter the jacobian.
You may ask that partial derivatives be stored as a sparse matrix by invoking
use_sparse_jacobians with a nonzero value. As with increasing dimension
jacobians tend to be sparsely occupied, doing so may eventually pay off in terms of
both memory consumption and speed.
octave:11> use_sparse_jacobians (1);
octave:12> [y, J] = D (@cumprod, reshape(1:9, 3, 3), 2)
y =
1 4 28
2 10 80
3 18 162
J =
Compressed Column Sparse (rows = 9, cols = 9, nnz = 18)
(1, 1) -> 1
(4, 1) -> 4
(7, 1) -> 28
(2, 2) -> 1
(5, 2) -> 5
(8, 2) -> 40
(3, 3) -> 1
(6, 3) -> 6
(9, 3) -> 54
(4, 4) -> 1
(7, 4) -> 7
(5, 5) -> 2
(8, 5) -> 16
(6, 6) -> 3
(9, 6) -> 27
(7, 7) -> 4
(8, 8) -> 10
(9, 9) -> 18
This is a best effort service, however, and there is no guarantee as to whether the returned jacobian will in fact be sparse. It certainly helps when the involved operands are:
octave:13> A = rand(6); b = rand(6); octave:14> x = sparse (A) \ gradinit (b); octave:15> spy (x.J, 0.5), issparse (x.J) ans = 1

Although primarily designing for functions with a real domain, the author does not think fit to impose any restriction here. Users should bear in mind though, when working with complex input, that what they get is the directional derivative along the real axis. It then may – or may not, for that matter – coincide with the derivative, depending on whether the function is locally holomorphic or not.
octave:16> [z, dz] = D (@abs, 1 + i)
z = 1.4142
dz = 0.70711
Beware that operator overloading is frail when it comes to interfacing with low-level routines. If you are in the habit of writing good portions of code in C++ or Fortran as DLD-functions – and there may well be good a reasons for it –, you will definitely run into trouble. The same caveat applies even to some functions of the Octave core API, in which case you should incur an error message like the one below:
octave:17> gamma (gradinit (4)) error: AD-rule unknown or function not overloaded
One might consider adding rules as the need arises. On the other hand,
balancing the benefit against the extra effort, it is often more reasonable to
fall back on numerical differentiation for less common operations and use
numgradient instead. In any event, the algebra provided by the
extension makes no claim for completeness and there certainly would be no point
in trying.