Numerical derivative with C++

Question

This question appeared a thousand times on different platforms. However I still need to understand something. Here is a complete example:

#include <iostream>
#include <iomanip>
#include <cmath>

// function to derivative
double f (double var, void * params){
  (void)(params); 
  return pow (var, 1.5);
}

// Naive central step derivative
double derivative1(double var, double f(double,void*), double h){ 
return (f(var+h,NULL) - f(var-h,NULL) )/2.0/h; 
}

// Richardson 5-point rule
double gderivative(double var, double f(double,void*), double h0){
 return(4.0*derivative1(var,f,h0) - derivative1(var,f,2.0*h0))/3.0;
}


int main (void){


  for (int i=10;i>0;i--){
  double h0=pow(10,-i);
  double x=2.0;
  double exact = 1.5 * sqrt(x);
  double test1=derivative1(x,f,h0);
  double gtest=gderivative(x,f,h0);
  std::cout << "h0 = " << h0 << std::endl;
  std::cout << "Exact      = "  << std::scientific<<std::setprecision(15) << exact << std::endl;
  std::cout << "Naive step = "  << std::setprecision(15) << test1 <<", diff = " << std::setprecision(15)<< exact-test1 << ", percent error = " << (exact-test1)/exact*100.0 << std::endl;
  std::cout << "Richardson = "  << std::setprecision(15) << gtest <<", diff = " << std::setprecision(15)<< exact-gtest << ", percent error = " << (exact-gtest)/exact*100.0 << std::endl;



}


 return 0;
}

The output is

h0 = 1e-10
Exact      = 2.121320343559643e+00
Naive step = 2.121318676273631e+00, diff = 1.667286011475255e-06, percent error = 7.859661632610832e-05
Richardson = 2.121318306199290e+00, diff = 2.037360352868944e-06, percent error = 9.604208808228318e-05
h0 = 1.000000000000000e-09
Exact      = 2.121320343559643e+00
Naive step = 2.121320674675076e+00, diff = -3.311154328500265e-07, percent error = -1.560893119491818e-05
Richardson = 2.121320748689944e+00, diff = -4.051303013064000e-07, percent error = -1.909802555452698e-05
h0 = 1.000000000000000e-08
Exact      = 2.121320343559643e+00
Naive step = 2.121320341608168e+00, diff = 1.951474537520426e-09, percent error = 9.199339191957163e-08
Richardson = 2.121320341608168e+00, diff = 1.951474537520426e-09, percent error = 9.199339191957163e-08
h0 = 1.000000000000000e-07
Exact      = 2.121320343559643e+00
Naive step = 2.121320341608168e+00, diff = 1.951474537520426e-09, percent error = 9.199339191957163e-08
Richardson = 2.121320340868019e+00, diff = 2.691623368633600e-09, percent error = 1.268843424240664e-07
h0 = 1.000000000000000e-06
Exact      = 2.121320343559643e+00
Naive step = 2.121320343606570e+00, diff = -4.692690680485612e-11, percent error = -2.212155601454860e-09
Richardson = 2.121320343643577e+00, diff = -8.393419292929138e-11, percent error = -3.956695799581460e-09
h0 = 1.000000000000000e-05
Exact      = 2.121320343559643e+00
Naive step = 2.121320343584365e+00, diff = -2.472244631235299e-11, percent error = -1.165427295665677e-09
Richardson = 2.121320343595468e+00, diff = -3.582467655860455e-11, percent error = -1.688791448560268e-09
h0 = 1.000000000000000e-04
Exact      = 2.121320343559643e+00
Naive step = 2.121320343340116e+00, diff = 2.195266191051815e-10, percent error = 1.034858406801534e-08
Richardson = 2.121320343561791e+00, diff = -2.148059508044753e-12, percent error = -1.012604963020456e-10
h0 = 1.000000000000000e-03
Exact      = 2.121320343559643e+00
Naive step = 2.121320321462283e+00, diff = 2.209735949776359e-08, percent error = 1.041679516479040e-06
Richardson = 2.121320343559311e+00, diff = 3.317346397579968e-13, percent error = 1.563812088849040e-11
h0 = 1.000000000000000e-02
Exact      = 2.121320343559643e+00
Naive step = 2.121318133840577e+00, diff = 2.209719065504601e-06, percent error = 1.041671557157002e-04
Richardson = 2.121320343601055e+00, diff = -4.141265108614789e-11, percent error = -1.952211093995174e-09
h0 = 1.000000000000000e-01
Exact      = 2.121320343559643e+00
Naive step = 2.121099269013200e+00, diff = 2.210745464426012e-04, percent error = 1.042155406248691e-02
Richardson = 2.121320759832334e+00, diff = -4.162726914280768e-07, percent error = -1.962328286210455e-05

I believe the standard GSL gsl_deriv_central employs the Richardson procedure.

Now the common agument that is given to choose h0 is that theoretically it should be chosen as small as possible to improve the precision of derivative however numerically it should not be too small so that we hit floating point round off thus ruining the precision. So often it is said that the optimal choice should be somewhat around 1e-6 - 1e-8. My question is :

What is the optimal choice for h0 in a generic derivative?

Should I have to check case by case? Often it might not be possible to have an exact result to check with. What should one do in that case? Now in this particular case it seems like the best choice for Naive step is h0 = 1.000000000000000e-05 whereas for Richardson h0 = 1.000000000000000e-03. This confuses me since these are not small. Any suggestion on any other good options (easy algorithm/library) which is efficient as well as precise(double)?

Answer 1

Personal opinion warning

In my experience i find it better to use a stepsize that is small compared to the variable it affects.

For example, I usually use something like this:

auto dfdx(std::function<double (double)> f, double x, double h0=1e-5){
    h0=std::abs(x)*h0+h0; // this is 1e-5 if x<<1e-5, and |x|*1e-5 otherwise
    
    return ( f(x+h0) - f(x-h0) )/2/h0;
}

This should work, since finite differences are motivated by Taylor expansion. That is, as long as x<<h0, the finite difference should be an good approximation.

Answer 2

This is quite as expected. The central difference has a second order error O(h^2) to the exact derivative. The function evaluation has an error of magnitude mu, the machine constant (for mildly scaled test examples). Thus the evaluation error of the central difference is of magnitude mu/h. The overall error is smallest if these two influences are about equal, thus h=mu^(1/3) gives h=1e-5, with an error of about 1e-10.

The same calculation for the Richardson extrapolation gives error order O(h^4) towards the exact derivative, resulting in h=mu^(1/5)=1e-3 as optimal step size, with an error of about 1e-12.

loglog plot of the errors of both methods over a larger sample of step sizes

In practical applications you would want to scale h so that the indicated sizes are relative to the size of x. The exact optimum depends also on the magnitudes of the derivatives of the function, if they grow too wildly.

To get more precise values for the derivatives, you would need a higher or multi-precision data type, or employ automatic/algorithmic differentiation, where the first and second derivatives can be evaluated with the same precision as the function values.