ComputeFiniteGradient/ComputeFiniteHessian: step size ignores `accuracy`, so higher-order stencils are no more accurate (and the Hessian is unusable)

[tesch1]

Summary

include/cppoptlib/utils/derivatives.h selects the finite-difference stencil from the accuracy argument but always uses h = sqrt(eps). The step is never matched to the order of the stencil that accuracy picked.

Two consequences:

1. ComputeFiniteGradient: accuracy 1, 2 and 3 cost 2x, 3x and 4x the function evaluations of accuracy 0 and return no extra accuracy. All four land at roughly 1e-9.
2. ComputeFiniteHessian: the same step is used in a second difference, so the divisor is h*h ≈ eps. The result is dominated by cancellation. On a smooth test function it returns -1.0 where the exact second derivative is -0.84147..., a 16% error. As a direct result, IsHessianCorrect returns false for a function whose analytic Hessian is exactly right.

Code

include/cppoptlib/utils/derivatives.h:68-73 (gradient):

  const int innerSteps = 2 * (accuracy + 1);
  for (index_t d = 0; d < x0.rows(); d++) {
    // Compute a coordinate-dependent step size.
    ScalarType h =
        std::sqrt(machine_eps) * std::max(std::abs(x0[d]), ScalarType(1));
    ScalarType ddVal = dd[accuracy] * h;

h does not read accuracy, while the stencils at lines 52-63 do:

  static const std::array<std::vector<ScalarType>, 4> coeff = {...};   // 2,4,6,8 points
  static const std::array<ScalarType, 4> dd = {2, 12, 60, 840};

Same sqrt(machine_eps) in ComputeFiniteHessian at lines 111, 123, 150 and 161, feeding divisors hi*hi (line 119), 4*hi*hj (line 141), hi*hi (line 157) and 600*h*h (line 246).

Why the step must follow the order

For a derivative of order d approximated to truncation order p, total error is roughly C1*h^p from truncation plus C2*eps/h^d from subtractive cancellation. Minimizing gives h* ~ eps^(1/(p+d)).

Observation 5.5.6 For computing with a finite-difference method of order m in the presence of roundoff, the optimal spacing of nodes satisfies
h_opt ≈ eps_mach^(1/(m+1)), (5.5.5)
and the optimum total error is roughly eps_mach^(m/(m+1)). [...] Higher-order finite-difference methods are both more efficient and less vulnerable to roundoff than low-order methods.

Driscoll & Braun, Fundamentals of Numerical Computation, §5.5 "Convergence of finite differences": https://tobydriscoll.net/fnc-julia/localapprox/fd-converge.html

Nocedal & Wright, Numerical Optimization (2nd ed.), §8.1, pp. 196-197 gives the two cases an optimization library cares about. For forward differences, eq. (8.6) is eps ≈ sqrt(u); for the central-difference formula (8.7),

"the same assumptions that were used to derive (8.6) lead to an optimal choice of ε of about u^{1/3} and an error of about u^{2/3}."

sqrt(eps) is the step for a forward difference. Every stencil in this file is a central difference, so even accuracy = 0 is using a step two orders of magnitude too small.

Measured

f = sin, x0 = 1, IEEE double, using this file's own coefficients, offsets and divisors:

accuracy	points	order p	err at h = sqrt(eps) = 1.49e-08	h* = eps^(1/(p+1))	err at h*
0	2	2	5.455e-10	6.055e-06	5.037e-12
1	4	4	6.963e-10	7.401e-04	9.237e-14
2	6	6	5.455e-10	5.805e-03	4.441e-16
3	8	8	6.874e-10	1.823e-02	1.776e-15

The h = sqrt(eps) column is flat, which is the point. It is round-off noise, so the exact values move a little with compiler and libm, but the magnitude does not. accuracy = 3 costs four times as much as accuracy = 0 and is no better. Stepped correctly it reaches 1.8e-15, close to six orders better than as shipped.

Hessian diagonal, (f(x+h) - 2f(x) + f(x-h)) / h^2, exact f''(1) = -0.8414709848078965:

h	result	error
sqrt(eps) = 1.490e-08	-1.0	1.585e-01
eps^(1/4) = 1.221e-04	-0.84147098660469055	1.797e-09

Against the library itself, f(x) = sin(x0)*cos(x1) at x = (1, 0) with analytic gradient and Hessian supplied:

=== gradient, |fd - exact|_inf ===
  accuracy=0 (2 evals/dim): 5.4551e-10
  accuracy=1 (4 evals/dim): 1.8626e-09
  accuracy=2 (6 evals/dim): 1.3659e-09
  accuracy=3 (8 evals/dim): 6.8743e-10
=== hessian, |fd - exact|_inf ===
  accuracy=0: 1.5853e-01   H(0,0)=-1 (exact -0.8414709848078965)
  accuracy=1: 1.5853e-01   H(0,0)=-1 (exact -0.8414709848078965)
=== validators on an ANALYTICALLY EXACT function ===
  IsGradientCorrect = true
  IsHessianCorrect  = FALSE

IsHessianCorrect uses tolerance = 1e-1 with scale >= 1 (lines 291, 301-307). The finite-difference error of 0.1585 exceeds that, so a correct user Hessian is reported as wrong.

Repro

Self-contained, no Eigen and no cppoptlib. c++ -std=c++17 -O2 standalone.cc -o fd && ./fd

// Replicates the stencils/divisors of ComputeFiniteGradient and the
// accuracy==0 diagonal of ComputeFiniteHessian.
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdio>
#include <limits>
#include <vector>

using S = double;
static const std::array<std::vector<S>, 4> coeff = {
    {{1, -1}, {1, -8, 8, -1}, {-1, 9, -45, 45, -9, 1},
     {3, -32, 168, -672, 672, -168, 32, -3}}};
static const std::array<std::vector<S>, 4> coeff2 = {
    {{1, -1}, {-2, -1, 1, 2}, {-3, -2, -1, 1, 2, 3},
     {-4, -3, -2, -1, 1, 2, 3, 4}}};
static const std::array<S, 4> dd = {2, 12, 60, 840};

S Grad(S x0, int accuracy, S h) {
  S g = 0, x = x0;
  for (int s = 0; s < 2 * (accuracy + 1); ++s) {
    S tmp = x; x += coeff2[accuracy][s] * h;
    g += coeff[accuracy][s] * std::sin(x); x = tmp;
  }
  return g / (dd[accuracy] * h);
}

int main() {
  constexpr S eps = std::numeric_limits<S>::epsilon();
  const S x0 = 1.0, exact = std::cos(x0);
  const S h_lib = std::sqrt(eps);  // what the library uses, for every order

  printf("acc pts order   h_lib      err(h_lib)   h*=eps^(1/(p+1))   err(h*)\n");
  for (int a = 0; a < 4; ++a) {
    const int p = 2 * (a + 1);
    const S h_opt = std::pow(eps, S(1) / S(p + 1));
    printf(" %d   %d    %d   %9.3e   %9.3e      %9.3e     %9.3e\n", a,
           2 * (a + 1), p, h_lib, std::fabs(Grad(x0, a, h_lib) - exact), h_opt,
           std::fabs(Grad(x0, a, h_opt) - exact));
  }

  auto D2 = [&](S h) {
    return (std::sin(x0 + h) - 2 * std::sin(x0) + std::sin(x0 - h)) / (h * h);
  };
  const S exact2 = -std::sin(x0);
  printf("\nHessian diagonal, f''(1) = %.17g\n", exact2);
  printf("  h = sqrt(eps) = %.3e -> %.17g  (err %.3e)\n",
         h_lib, D2(h_lib), std::fabs(D2(h_lib) - exact2));
  const S h4 = std::pow(eps, S(1) / S(4));
  printf("  h = eps^(1/4) = %.3e -> %.17g  (err %.3e)\n",
         h4, D2(h4), std::fabs(D2(h4) - exact2));
  return 0;
}

Suggested fix

Gradient, derivatives.h:68-73. The stencil order equals innerSteps, so h ~ eps^(1/(order+1)). Hoist it out of the loop since it does not depend on d:

  const int innerSteps = 2 * (accuracy + 1);
  const int order = 2 * (accuracy + 1);
  const ScalarType h_rel =
      std::pow(machine_eps, ScalarType(1) / ScalarType(order + 1));
  for (index_t d = 0; d < x0.rows(); d++) {
    ScalarType h = h_rel * std::max(std::abs(x0[d]), ScalarType(1));
    ScalarType ddVal = dd[accuracy] * h;

Hessian: a second derivative has round-off ~eps/h^2, so the exponent is 1/(p+2), not 1/(p+1).

• accuracy == 0, both the diagonal (line 119) and the off-diagonal (line 141) are 2nd order, verified by a convergence sweep. Use eps^(1/4) ≈ 1.22e-4 at lines 111 and 123.
• accuracy != 0, the 600*h*h mixed formula (line 246) is 4th order, also verified by sweep, so it wants eps^(1/6) ≈ 2.6e-3 at lines 161-162.
• Note the diagonal in the accuracy != 0 branch (lines 149-157) is the same 2nd-order 3-point formula as the accuracy == 0 branch. accuracy > 0 never improves the diagonal at all. That is arguably a separate bug, but at minimum its hi should use the 2nd-order exponent eps^(1/4), not the 4th-order one.

Minor: the file uses std::sqrt/std::abs without including <cmath> and currently gets it transitively. std::pow should come with an explicit #include <cmath>.

A regression test is easy to add. IsHessianCorrect on any smooth function with an exact analytic Hessian currently fails and should pass.