Improve numeric differentation near zero.
Before this change, the default step size
for a function F(x) at x was
step_size = |x| * relative_step_size
if step_size was exactly zero, then to prevent
division by zero we would fall back to relative_step_size.
This however is not good enough, as values of x say 1e-64
would lead to step sizes ~ 1e-70 and dividing by such numbers
leads to inaccurate results. For even smaller numbers, like
1e-300, which I have observed can occur as the optimization
algorithm makes progress, this leads to NaNs.
The key change in this CL is to change the fallback mechanism
to be
step_size = max(|x| * relative_step_size, min_step_size)
where
min_step_size = sqrt(DBL_EPSILON)
This is the recommended minimum value for the step size
for double precision arithmetic on the interwebs.
This results in a small loss of precision in the transcendental
functions test, but that is unavoidable as we are not taking
sufficiently small steps anymore.
On the whole though this will improve the numerical performance
of the algorithm.
To validate this approach, one of the parameter values for the
EasyFunctorTest has been set to 1e-64, which causes the test
to start failing without the corrected fallback logic.
This change should also address some if not all of
https://github.com/ceres-solver/ceres-solver/issues/121
Change-Id: I4a9013ef358626c1ba7b8abad60b3904163d63f6
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);
