I've been working recently on the AMPL interface for the GNU Scientific Library. Adding a function to AMPL normally requires providing first and second partial derivatives for it. Computing the derivatives is straightforward and there are plenty of tools, such as Wolfram Alpha, that automate this. However, transferring this to code is a tedious and error-prone process.

For example, consider the function gsl_sf_bessel_I0_scaled(double x) which is defined as \(e^{-|x|} I_0(x)\), where \(I_0(x)\) is the regular modified cylindrical Bessel function of zeroth order. It has a simple derivative:

\( e^{-|x|} I_1(x)- \frac{x}{|x|} e^{-|x|} I_0(x) \)

The derivative can be expressed even in a more compact form using \(I_0(x)\) and \(I_1(x)\) scaled by \(e^{-|x|}\), gsl_sf_bessel_I0_scaled and gsl_sf_bessel_I1_scaled.

Transferring this into C results in the following small program:

#include <math.h>
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

double I0_scaled_deriv(double x) {
  return gsl_sf_bessel_I1_scaled(x) -
    x * gsl_sf_bessel_I0_scaled(x) / abs(x);
}

int main() {
  printf("%g\n", I0_scaled_deriv(1.5));
  return 0;
}

Let's compile and run the program:

$ gcc -pedantic -Wall -Wextra test.c -lgsl -lgslcblas -otest
$ ./test 
-0.332111

It is easy to check that the printed value is incorrect, -0.148394 should be printed instead. Can you spot the error?

Show/Hide the answer

This small example demonstrates how important it is to test the code. Unfortunately it is rarely discussed in the context of OR. As Tim Hopper put it in his recent blog post:

Code testing was never mentioned in any of my classes in college, and I never hear operations researchers talk about it. Again, I have no doubt that the code behind much published work is full of mistakes. Operations researchers need good testing practices?

Returning to the question of derivatives, how do you actually test that the implementation is correct? One possibility is to use numerical differentiation. The simplest method is directly based on the definition of the derivative:

\( f'(x) \approx \frac{f(x + h) - f(x)}{h} \) for some small \(h\)

Although trivial to implement, this method suffers from several issues such as the difficulty of selecting \(h\) especially for \(x = 0\). So I used the Ridders' method described in the book Numerical Recipes in C: The Art of Scientific Computing. This methods starts from large \(h\) and successfully decreases it estimating error at each iteration.

Numerical differentiation proved to be very useful for testing. Although it took some time to implement the test infrastructure, it paid off by catching a number of bugs and corner cases. And now that the test infrastructure is in place, adding tests for new functions requires minimal amount of code.

Last modified on 2012-07-17