Numerical differentiation by finite differences

Numerical differentiation is known to be ill-conditioned unless using a Chebyshev series, but this requires global information about the function and a priori knowledge of a compact domain on which the function will be evaluated.

For this reason, simple finite differences are often useful. The unit roundoff gives a natural choice of step-see here and here for more details. I felt that this technique is better than the technique used in the GSL which requires specifying the stepsize, so I felt it should be given life in Boost.Math, esp. since spectral element methods are getting so popular.

Regrettably, the pull-request was panned because the implementation is awkward, requiring a bunch of if checks on the order of accuracy and quite a bit of code duplication. Can it be cleaned up?

#ifndef BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP
#define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP
#include <boost/math/constants/constants.hpp>

namespace boost { namespace math { namespace tools {
template<class F, class Real, size_t order=6>
Real finite_difference_derivative(const F f, Real x, Real* error = nullptr)
{
static_assert(order == 1 || order == 2 || order == 4 || order == 6 || order == 8,
"Order of accuracy must be one of 1, 2, 4, 6, or 8.\n");
// These using declarations allow for ADL so the boost.multiprecision types work.
using std::sqrt;
using std::pow;
using std::abs;
using std::max;
using std::nextafter;
using boost::math::constants::half;

const Real eps = std::numeric_limits<Real>::epsilon();
// TODO: static if in C++17!
// Error bound ~eps^1/2
if (order == 1)
{
// Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|).
// Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1.
// This approximation will get better as we move to higher orders of accuracy.
Real h = 2*sqrt(eps);

// Redefine h so that x + h is representable. Not using this trick leads to large error.
// The compiler flag -ffast-math evaporates these operations . . .
Real temp = x + h;
h = temp - x;
// Handle the case x + h == x:
if (h == 0)
{
h = nextafter(x, std::numeric_limits<Real>::max()) - x;
}
Real yh = f(x+h);
Real y0 = f(x);
Real diff = yh - y0;
if (error)
{
Real ym = f(x-h);
Real ypph = abs(yh - 2*y0 + ym)/h;
// h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h,
*error = ypph*half<Real>() + (abs(yh) + abs(y0))*eps/h;
}
return diff/h;
}

// Error bound ~eps^2/3
if (order == 2)
{
// See the previous discussion to understand determination of h and the error bound.
// Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}]
Real h = pow(3*eps, boost::math::constants::third<Real>());
Real temp = x + h;
h = temp - x;
if (h == 0)
{
h = nextafter(x, std::numeric_limits<Real>::max()) - x;
}
Real yh = f(x+h);
Real ymh = f(x-h);
Real diff = yh - ymh;
if (error)
{
Real yth = f(x+2*h);
Real ymth = f(x-2*h);
*error = eps*(abs(yh) + abs(ymh))/(2*h) + abs((yth - ymth)*half<Real>() - diff)/(6*h);
}

return diff/(2*h);
}

// Error bound ~eps^4/5
if (order == 4)
{
Real h = pow(11.25*eps, (Real) 1/ (Real) 5);
Real temp = x + h;
h = temp - x;
if (h == 0)
{
h = nextafter(x, std::numeric_limits<Real>::max()) - x;
}
Real ymth = f(x-2*h);
Real yth = f(x+2*h);
Real yh = f(x+h);
Real ymh = f(x-h);
Real y2 = ymth - yth;
Real y1 = yh - ymh;
if (error)
{
// Mathematica code to extrace the remainder:
// Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}]
Real y_three_h = f(x+3*h);
Real y_m_three_h = f(x-3*h);
// Error from fifth derivative:
*error = abs(half<Real>()*(y_three_h - y_m_three_h) + 2*(ymth - yth) + 5*half<Real>()*(yh - ymh) )/(30*h);
// Error from function evaluation
*error += eps*(abs(yth) + abs(ymth) + 8*(abs(ymh) + abs(yh)))/(12*h);
}
return (y2+8*y1)/(12*h);
}

// Error bound ~eps^6/7
if (order == 6)
{
// Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h
Real h = pow(eps/168, (Real) 1/ (Real) 7);
Real temp = x + h;
h = temp - x;
if (h == 0)
{
h = nextafter(x, std::numeric_limits<Real>::max()) - x;
}

Real yh = f(x+h);
Real ymh = f(x-h);
Real y1 = yh - ymh;
Real y2 = f(x-2*h) - f(x + 2*h);
Real y3 = f(x+3*h) - f(x - 3*h);

if (error)
{
// Mathematica code to generate fd scheme for 7th derivative:
// Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}]
// Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
// Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}]
Real y7 = half<Real>()*(f(x+4*h) - f(x-4*h) - 6*y3 - 14*y1 - 14*y2);
*error = abs(y7)/(140*h) + 5*(abs(yh) + abs(ymh))*eps/h;
}
return (y3 + 9*y2 + 45*y1)/(60*h);
}

// Error bound ~eps^8/9.
// In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations.
if (order == 8)
{
// Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions.
// Mathematica code to get the error:
// Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}]
// If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
Real h = pow(551.25*eps, (Real)1 / (Real) 9);
Real temp = x + h;
h = temp - x;
if (h == 0)
{
h = nextafter(x, std::numeric_limits<Real>::max()) - x;
}

Real yh = f(x+h);
Real ymh = f(x-h);
Real y1 = yh - ymh;
Real y2 = f(x-2*h) - f(x + 2*h);
Real y3 = f(x+3*h) - f(x - 3*h);
Real y4 = f(x-4*h) - f(x + 4*h);

Real tmp1 = 3*y4/8 + 4*y3;
Real tmp2 = 21*y2 + 84*y1;

if (error)
{
// Mathematica code to generate fd scheme for 7th derivative:
// Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}]
// Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
// Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h]  - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}]
Real f9 = (f(x+5*h) - f(x-5*h))*half<Real>() + 4*y4 + 27*half<Real>()*y3 + 24*y2 + 21*y1;
*error = abs(f9)/(630*h) + 7*(abs(yh)+abs(ymh))*eps/h;
}
return (tmp1 + tmp2)/(105*h);
}
}

}}}
#endif

Can you make a generic implementation that handles any order?

You are handling each of the supported orders completely separately. That's not very generic, and then you might wonder if it is a good idea to make the order a template argument in the first place. But there is a structure to how you calculate an n'th order derivative. So try to write generic code that handles any given order, or at least so that there are only two versions, one for odd and one for even orders.

In particular, this would then also get rid of all the magic constants that you have in the code right now, as it would force you to write an expression to derive those constants.

Guarding against compiler optimizations

You already mentioned that -ffast-math will optimize away the following:

Real temp = x + h;
h = temp - x;

But what you actually want to ensure is that x + h != x. But that also might be optimized away. So how about checking that x + h >= nextafter(x, std::numeric_limits<Real>::max)? And you can avoid an if-statement by writing:

h = std::max(x + h, std::nextafter(x, std::numeric_limit<Real>::max)) - x;

Fractions

I'm not sure what good std::boost::constants::half<Real> and related constants bring here, when it obviously is limited in the number of fractions that are pre-defined, and that you have to fall back to writing fractions the normal way for anyway. Try to be consistent, and since you can't be consistent using std::boost::constants, I would not use it at all. Instead of writing (Real)1 / (Real)/5, I would write Real(1) / 5.

Also, just divide Real values by two by writing / 2 instead of * half<Real>(), the latter is very hard to read and brings no benefit.