NumericalDiff: Divide by the step actually applied to x

libeigen/eigen!2987

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

unsupported/Eigen/src/NumericalDiff/NumericalDiff.h

diff --git a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
index f07d3cc..1944f98 100644
--- a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
+++ b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
@@ -19,6 +19,21 @@
 
 namespace Eigen {
 
+namespace internal {
+
+// Keeps -ffast-math from folding the rounded evaluation points back into x and h, e.g. (x + h) - x
+// into h. EIGEN_OPTIMIZATION_BARRIER is an asm operand constraint restricted to plain types (see
+// Macros.h), so class-type scalars go unguarded.
+template <typename Scalar, std::enable_if_t<std::is_floating_point<Scalar>::value, int> = 0>
+EIGEN_STRONG_INLINE void numerical_diff_barrier(Scalar& x) {
+  EIGEN_UNUSED_VARIABLE(x);
+  EIGEN_OPTIMIZATION_BARRIER(x)
+}
+template <typename Scalar, std::enable_if_t<!std::is_floating_point<Scalar>::value, int> = 0>
+EIGEN_STRONG_INLINE void numerical_diff_barrier(Scalar&) {}
+
+}  // namespace internal
+
 enum NumericalDiffMode { Forward, Central };
 
 /**
@@ -55,7 +70,11 @@
   enum { InputsAtCompileTime = Functor::InputsAtCompileTime, ValuesAtCompileTime = Functor::ValuesAtCompileTime };
 
   /**
-   * return the number of evaluations of the functor
+   * Computes the Jacobian of the functor at \a _x into \a jac and returns the number of functor evaluations.
+   *
+   * The step along coordinate \c j is <tt>h = eps * max(|x[j]|, 1)</tt> with <tt>eps = sqrt(max(epsfcn, epsilon))</tt>
+   * and \c epsilon the machine precision NumTraits<Scalar>::epsilon(); the difference quotient divides by the
+   * representable step <tt>fl(x[j] + h) - x[j]</tt> actually applied.
    */
   int df(const InputType& _x, JacobianType& jac) const {
     using std::abs;
@@ -89,24 +108,35 @@
     for (int j = 0; j < n; ++j) {
       const Scalar x_abs = abs(x[j]);
       h = numext::maxi(x_abs, Scalar(1)) * eps;
+      // The functor is evaluated at fl(x[j] + h), so divide by that representable step: the rounding
+      // of x[j] + h perturbs h by up to ulp(x[j]) <= epsilon/eps * h <= sqrt(epsilon) * h, comparable
+      // to the error of the difference quotient itself.
+      Scalar x_plus = _x[j] + h;
+      internal::numerical_diff_barrier(x_plus);
+      h = x_plus - _x[j];
       switch (mode) {
         case Forward:
-          x[j] += h;
+          x[j] = x_plus;
           Functor::operator()(x, val2);
           nfev++;
           x[j] = _x[j];
           jac.col(j) = (val2 - val1) / h;
           break;
-        case Central:
-          x[j] += h;
+        case Central: {
+          x[j] = x_plus;
           Functor::operator()(x, val2);
           nfev++;
-          x[j] -= 2 * h;
+          // x[j] - h can round (a tie when x[j] < 0 has |x[j]| within h above a power of two), so
+          // divide by the separation of the two evaluation points rather than by 2*h.
+          Scalar x_minus = _x[j] - h;
+          internal::numerical_diff_barrier(x_minus);
+          x[j] = x_minus;
           Functor::operator()(x, val1);
           nfev++;
           x[j] = _x[j];
-          jac.col(j) = (val2 - val1) / (2 * h);
+          jac.col(j) = (val2 - val1) / (x_plus - x_minus);
           break;
+        }
         default:
           eigen_assert(false);
       }

unsupported/test/NumericalDiff.cpp

diff --git a/unsupported/test/NumericalDiff.cpp b/unsupported/test/NumericalDiff.cpp
index bfde6b6..d2ebd6a 100644
--- a/unsupported/test/NumericalDiff.cpp
+++ b/unsupported/test/NumericalDiff.cpp
@@ -59,6 +59,17 @@
   }
 };
 
+// Scaling by a power of two is exact, so the difference quotients of this map equal the scales
+// iff each divisor is the representable step actually applied to x.
+struct pow2_scaling_functor : Functor<double> {
+  pow2_scaling_functor() : Functor<double>(7, 7), scale(7) { scale << 0.25, 1.0, 8.0, -2.0, 0.5, 4.0, -0.125; }
+  int operator()(const VectorXd &x, VectorXd &fvec) const {
+    fvec = scale.cwiseProduct(x);
+    return 0;
+  }
+  VectorXd scale;
+};
+
 void test_forward() {
   VectorXd x(3);
   MatrixXd jac(15, 3);
@@ -115,8 +126,29 @@
   VERIFY_IS_APPROX(jac, actual_jac);
 }
 
+template <NumericalDiffMode Mode>
+void test_step_alignment(double epsfcn) {
+  pow2_scaling_functor functor;
+  NumericalDiff<pow2_scaling_functor, Mode> numDiff(functor, epsfcn);
+
+  // Non-dyadic coordinates, so that x + h rounds; the first two take the absolute step eps, and the
+  // last two are negative with |x| at or just above a power of two, where x - h is a rounding tie.
+  VectorXd x(7);
+  x << 0.0, 0.082, 1.13, -2350.7, 7.7e8 + 0.3, -1.0 - std::ldexp(1.0, -27), -1.0;
+  MatrixXd jac(7, 7);
+  numDiff.df(x, jac);
+
+  MatrixXd expected = functor.scale.asDiagonal();
+  VERIFY_IS_EQUAL(jac, expected);
+}
+
 EIGEN_DECLARE_TEST(NumericalDiff) {
   CALL_SUBTEST(test_forward());
   CALL_SUBTEST(test_central());
   CALL_SUBTEST(test_small());
+  CALL_SUBTEST(test_step_alignment<Forward>(0.0));
+  CALL_SUBTEST(test_step_alignment<Central>(0.0));
+  // A non-dyadic eps also misaligns the absolute step taken for |x| < 1.
+  CALL_SUBTEST(test_step_alignment<Forward>(1e-10));
+  CALL_SUBTEST(test_step_alignment<Central>(1e-10));
 }