Complex Step

Finite differences are cool, but you are limited by subtractive cancellation (section 3.4: subtractive cancellation exercises). What if you want to make a numerical estimate of the derivative of a function, like a 2D Gaussian:


A finite difference approximation would be (in Octave):

dfdx = ( f(x+h) - f(x) )/ h;

Complex step is even cooler, because you don't have any subtraction (there's no difference), so you can choose a very small step size without loosing accuracy due to subtractive cancelation:

complex_step = complex( 0, 1e-320 );

The derivative is approximated by just the imaginary part:

x = x + complex_step;
dfdx = imag( f(x) )/imag( complex_step );

So the derivative with respect to x looks like this:

And the estimated derivative with respect to y:


This approach is really useful for design sensitivity analysis, and since modern Fortran supports complex types we can even use this method for serious number crunching!