diff --git a/include/ceres/internal/numeric_diff.h b/include/ceres/internal/numeric_diff.h
index d220c7c..da04a8d 100644
--- a/include/ceres/internal/numeric_diff.h
+++ b/include/ceres/internal/numeric_diff.h
@@ -133,23 +133,18 @@
     ParameterVector x(x_plus_delta);
     ParameterVector step_size = x.array().abs() * relative_step_size;
 
-    // To handle cases where a parameter is exactly zero, instead use
-    // the mean step_size for the other dimensions. If all the
-    // parameters are zero, there's no good answer. Take
-    // relative_step_size as a guess and hope for the best.
-    const double fallback_step_size =
-        (step_size.sum() == 0)
-        ? relative_step_size
-        : step_size.sum() / step_size.rows();
+    // It is not a good idea to make the step size arbitrarily
+    // small. This will lead to problems with round off and numerical
+    // instability when dividing by the step size. The general
+    // recommendation is to not go down below sqrt(epsilon).
+    const double min_step_size =
+        std::sqrt(std::numeric_limits<double>::epsilon());
 
     // For each parameter in the parameter block, use finite differences to
     // compute the derivative for that parameter.
-
     ResidualVector residuals(num_residuals_internal);
     for (int j = 0; j < parameter_block_size_internal; ++j) {
-      const double delta =
-          (step_size(j) == 0.0) ? fallback_step_size : step_size(j);
-
+      const double delta = std::max(min_step_size, step_size(j));
       x_plus_delta(j) = x(j) + delta;
 
       if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
diff --git a/internal/ceres/numeric_diff_test_utils.cc b/internal/ceres/numeric_diff_test_utils.cc
index 95d1da5..5595449 100644
--- a/internal/ceres/numeric_diff_test_utils.cc
+++ b/internal/ceres/numeric_diff_test_utils.cc
@@ -57,7 +57,9 @@
 void EasyFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
     const CostFunction& cost_function,
     NumericDiffMethod method) const {
-  double x1[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+  // The x1[0] is made deliberately small to test the performance near
+  // zero.
+  double x1[] = { 1e-64, 2.0, 3.0, 4.0, 5.0 };
   double x2[] = { 9.0, 9.0, 5.0, 5.0, 1.0 };
   double *parameters[] = { &x1[0], &x2[0] };
 
@@ -71,9 +73,12 @@
                                      &residuals[0],
                                      &jacobians[0]));
 
-  EXPECT_EQ(residuals[0], 67);
-  EXPECT_EQ(residuals[1], 4489);
-  EXPECT_EQ(residuals[2], 213);
+  double expected_residuals[3];
+  EasyFunctor functor;
+  functor(x1, x2, expected_residuals);
+  EXPECT_EQ(expected_residuals[0], residuals[0]);
+  EXPECT_EQ(expected_residuals[1], residuals[1]);
+  EXPECT_EQ(expected_residuals[2], residuals[2]);
 
   const double tolerance = (method == CENTRAL)? 3e-9 : 2e-5;
 
@@ -145,7 +150,7 @@
       x1x2 += x1[i] * x2[i];
     }
 
-    const double tolerance = (method == CENTRAL)? 3e-9 : 2e-5;
+    const double tolerance = (method == CENTRAL)? 2e-7 : 2e-5;
 
     for (int i = 0; i < 5; ++i) {
       ExpectClose( x2[i] * cos(x1x2),              dydx1[5 * 0 + i], tolerance);