Automatic Differentiation (AD) in Octave


Contents


Concept

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.

Octave-AD-Extension

License Information and Disclaimer

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.

Prerequisites

The maintainer will do his best to retain upward compatibility, though at the current pace of releases this seems an audacious promise.

New Features

Download and Installation

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.

Testsuite

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.

Using AD in Octave

User Interface and Class Gradient

The function 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;

With 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.

Sparse Storage Mode

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
  
  

Complex-valued Domains

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

Limitations

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.


Last Modified: Mon Feb 18 10:19:12 CET 2008