Adaptive numeric differentiation using Ridders' method.
This method numerically computes function derivatives in different
scales, extrapolating between intermediate results to conserve function
evaluations. Adaptive differentiation is essential to produce accurate
results for functions with noisy derivatives.
Full changelist:
-Created a new type of NumericDiffMethod (RIDDERS).
-Implemented EvaluateRiddersJacobianColumn in NumericDiff.
-Created unit tests with f(x) = x^2 + [random noise] and
f(x) = exp(x).
Change-Id: I2d6e924d7ff686650272f29a8c981351e6f72091
diff --git a/internal/ceres/gradient_checking_cost_function.cc b/internal/ceres/gradient_checking_cost_function.cc
index e3ad3d4..580fd26 100644
--- a/internal/ceres/gradient_checking_cost_function.cc
+++ b/internal/ceres/gradient_checking_cost_function.cc
@@ -86,7 +86,7 @@
class GradientCheckingCostFunction : public CostFunction {
public:
GradientCheckingCostFunction(const CostFunction* function,
- double relative_step_size,
+ const NumericDiffOptions& options,
double relative_precision,
const string& extra_info)
: function_(function),
@@ -97,7 +97,7 @@
new DynamicNumericDiffCostFunction<CostFunction, CENTRAL>(
function,
DO_NOT_TAKE_OWNERSHIP,
- relative_step_size);
+ options);
const vector<int32>& parameter_block_sizes =
function->parameter_block_sizes();
@@ -235,8 +235,11 @@
double relative_step_size,
double relative_precision,
const string& extra_info) {
+ NumericDiffOptions numeric_diff_options;
+ numeric_diff_options.relative_step_size = relative_step_size;
+
return new GradientCheckingCostFunction(cost_function,
- relative_step_size,
+ numeric_diff_options,
relative_precision,
extra_info);
}
diff --git a/internal/ceres/numeric_diff_cost_function_test.cc b/internal/ceres/numeric_diff_cost_function_test.cc
index cb74a86..13ab106 100644
--- a/internal/ceres/numeric_diff_cost_function_test.cc
+++ b/internal/ceres/numeric_diff_cost_function_test.cc
@@ -27,6 +27,7 @@
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
+// tbennun@gmail.com (Tal Ben-Nun)
#include "ceres/numeric_diff_cost_function.h"
@@ -71,6 +72,19 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
}
+TEST(NumericDiffCostFunction, EasyCaseFunctorRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<EasyFunctor,
+ RIDDERS,
+ 3, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new EasyFunctor));
+ EasyFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
TEST(NumericDiffCostFunction, EasyCaseCostFunctionCentralDifferences) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
@@ -97,7 +111,21 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
}
-TEST(NumericDiffCostFunction, TranscendentalCaseFunctorCentralDifferences) {
+TEST(NumericDiffCostFunction, EasyCaseCostFunctionRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<EasyCostFunction,
+ RIDDERS,
+ 3, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new EasyCostFunction, TAKE_OWNERSHIP));
+ EasyFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseFunctorCentralDifferences) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
new NumericDiffCostFunction<TranscendentalFunctor,
@@ -110,7 +138,8 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
}
-TEST(NumericDiffCostFunction, TranscendentalCaseFunctorForwardDifferences) {
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseFunctorForwardDifferences) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
new NumericDiffCostFunction<TranscendentalFunctor,
@@ -123,7 +152,28 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
}
-TEST(NumericDiffCostFunction, TranscendentalCaseCostFunctionCentralDifferences) {
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseFunctorRidders) {
+ NumericDiffOptions options;
+
+ // Using a smaller initial step size to overcome oscillatory function
+ // behavior.
+ options.ridders_relative_initial_step_size = 1e-3;
+
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<TranscendentalFunctor,
+ RIDDERS,
+ 2, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new TranscendentalFunctor, TAKE_OWNERSHIP, 2, options));
+ TranscendentalFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseCostFunctionCentralDifferences) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
new NumericDiffCostFunction<TranscendentalCostFunction,
@@ -136,7 +186,8 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
}
-TEST(NumericDiffCostFunction, TranscendentalCaseCostFunctionForwardDifferences) {
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseCostFunctionForwardDifferences) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
new NumericDiffCostFunction<TranscendentalCostFunction,
@@ -149,6 +200,26 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD);
}
+TEST(NumericDiffCostFunction,
+ TranscendentalCaseCostFunctionRidders) {
+ NumericDiffOptions options;
+
+ // Using a smaller initial step size to overcome oscillatory function
+ // behavior.
+ options.ridders_relative_initial_step_size = 1e-3;
+
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<TranscendentalCostFunction,
+ RIDDERS,
+ 2, /* number of residuals */
+ 5, /* size of x1 */
+ 5 /* size of x2 */>(
+ new TranscendentalCostFunction, TAKE_OWNERSHIP, 2, options));
+ TranscendentalFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS);
+}
+
template<int num_rows, int num_cols>
class SizeTestingCostFunction : public SizedCostFunction<num_rows, num_cols> {
public:
@@ -208,7 +279,8 @@
new EasyFunctor, TAKE_OWNERSHIP, 2));
}
-TEST(NumericDiffCostFunction, EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) {
+TEST(NumericDiffCostFunction,
+ EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) {
internal::scoped_ptr<CostFunction> cost_function;
cost_function.reset(
new NumericDiffCostFunction<EasyFunctor,
@@ -221,5 +293,66 @@
functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, CENTRAL);
}
+TEST(NumericDiffCostFunction, ExponentialFunctorRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<ExponentialFunctor,
+ RIDDERS,
+ 1, /* number of residuals */
+ 1 /* size of x1 */>(
+ new ExponentialFunctor));
+ ExponentialFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, ExponentialCostFunctionRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ cost_function.reset(
+ new NumericDiffCostFunction<ExponentialCostFunction,
+ RIDDERS,
+ 1, /* number of residuals */
+ 1 /* size of x1 */>(
+ new ExponentialCostFunction));
+ ExponentialFunctor functor;
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, RandomizedFunctorRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ NumericDiffOptions options;
+ // Larger initial step size is chosen to produce robust results in the
+ // presence of random noise.
+ options.ridders_relative_initial_step_size = 10.0;
+
+ cost_function.reset(
+ new NumericDiffCostFunction<RandomizedFunctor,
+ RIDDERS,
+ 1, /* number of residuals */
+ 1 /* size of x1 */>(
+ new RandomizedFunctor(kNoiseFactor, kRandomSeed), TAKE_OWNERSHIP,
+ 1, options));
+ RandomizedFunctor functor (kNoiseFactor, kRandomSeed);
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+TEST(NumericDiffCostFunction, RandomizedCostFunctionRidders) {
+ internal::scoped_ptr<CostFunction> cost_function;
+ NumericDiffOptions options;
+ // Larger initial step size is chosen to produce robust results in the
+ // presence of random noise.
+ options.ridders_relative_initial_step_size = 10.0;
+
+ cost_function.reset(
+ new NumericDiffCostFunction<RandomizedCostFunction,
+ RIDDERS,
+ 1, /* number of residuals */
+ 1 /* size of x1 */>(
+ new RandomizedCostFunction(kNoiseFactor, kRandomSeed),
+ TAKE_OWNERSHIP, 1, options));
+ RandomizedFunctor functor (kNoiseFactor, kRandomSeed);
+ functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function);
+}
+
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/numeric_diff_test_utils.cc b/internal/ceres/numeric_diff_test_utils.cc
index 5595449..bab186c 100644
--- a/internal/ceres/numeric_diff_test_utils.cc
+++ b/internal/ceres/numeric_diff_test_utils.cc
@@ -27,6 +27,7 @@
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
+// tbennun@gmail.com (Tal Ben-Nun)
#include "ceres/numeric_diff_test_utils.h"
@@ -56,7 +57,7 @@
void EasyFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
const CostFunction& cost_function,
- NumericDiffMethod method) const {
+ NumericDiffMethodType method) const {
// 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 };
@@ -80,7 +81,21 @@
EXPECT_EQ(expected_residuals[1], residuals[1]);
EXPECT_EQ(expected_residuals[2], residuals[2]);
- const double tolerance = (method == CENTRAL)? 3e-9 : 2e-5;
+ double tolerance = 0.0;
+ switch (method) {
+ default:
+ case CENTRAL:
+ tolerance = 3e-9;
+ break;
+
+ case FORWARD:
+ tolerance = 2e-5;
+ break;
+
+ case RIDDERS:
+ tolerance = 1e-13;
+ break;
+ }
for (int i = 0; i < 5; ++i) {
ExpectClose(x2[i], dydx1[5 * 0 + i], tolerance); // y1
@@ -106,7 +121,7 @@
void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
const CostFunction& cost_function,
- NumericDiffMethod method) const {
+ NumericDiffMethodType method) const {
struct {
double x1[5];
double x2[5];
@@ -150,7 +165,21 @@
x1x2 += x1[i] * x2[i];
}
- const double tolerance = (method == CENTRAL)? 2e-7 : 2e-5;
+ double tolerance = 0.0;
+ switch (method) {
+ default:
+ case CENTRAL:
+ tolerance = 2e-7;
+ break;
+
+ case FORWARD:
+ tolerance = 2e-5;
+ break;
+
+ case RIDDERS:
+ tolerance = 3e-12;
+ break;
+ }
for (int i = 0; i < 5; ++i) {
ExpectClose( x2[i] * cos(x1x2), dydx1[5 * 0 + i], tolerance);
@@ -161,5 +190,81 @@
}
}
+bool ExponentialFunctor::operator()(const double* x1,
+ double* residuals) const {
+ residuals[0] = exp(x1[0]);
+ return true;
+}
+
+void ExponentialFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
+ const CostFunction& cost_function) const {
+ // Evaluating the functor at specific points for testing.
+ double kTests[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+
+ // Minimal tolerance w.r.t. the cost function and the tests.
+ const double kTolerance = 2e-14;
+
+ for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
+ double *parameters[] = { &kTests[k] };
+ double dydx;
+ double *jacobians[1] = { &dydx };
+ double residual;
+
+ ASSERT_TRUE(cost_function.Evaluate(¶meters[0],
+ &residual,
+ &jacobians[0]));
+
+
+ double expected_result = exp(kTests[k]);
+
+ // Expect residual to be close to exp(x).
+ ExpectClose(residual, expected_result, kTolerance);
+
+ // Check evaluated differences. dydx should also be close to exp(x).
+ ExpectClose(dydx, expected_result, kTolerance);
+ }
+}
+
+bool RandomizedFunctor::operator()(const double* x1,
+ double* residuals) const {
+ double random_value = static_cast<double>(rand()) /
+ static_cast<double>(RAND_MAX);
+
+ // Normalize noise to [-factor, factor].
+ random_value *= 2.0;
+ random_value -= 1.0;
+ random_value *= noise_factor_;
+
+ residuals[0] = x1[0] * x1[0] + random_value;
+ return true;
+}
+
+void RandomizedFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect(
+ const CostFunction& cost_function) const {
+ double kTests[] = { 0.0, 1.0, 3.0, 4.0, 50.0 };
+
+ const double kTolerance = 2e-4;
+
+ // Initialize random number generator with given seed.
+ srand(random_seed_);
+
+ for (int k = 0; k < CERES_ARRAYSIZE(kTests); ++k) {
+ double *parameters[] = { &kTests[k] };
+ double dydx;
+ double *jacobians[1] = { &dydx };
+ double residual;
+
+ ASSERT_TRUE(cost_function.Evaluate(¶meters[0],
+ &residual,
+ &jacobians[0]));
+
+ // Expect residual to be close to x^2 w.r.t. noise factor.
+ ExpectClose(residual, kTests[k] * kTests[k], noise_factor_);
+
+ // Check evaluated differences. (dy/dx = ~2x)
+ ExpectClose(dydx, 2 * kTests[k], kTolerance);
+ }
+}
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/numeric_diff_test_utils.h b/internal/ceres/numeric_diff_test_utils.h
index 3ae83fa..2a551d3 100644
--- a/internal/ceres/numeric_diff_test_utils.h
+++ b/internal/ceres/numeric_diff_test_utils.h
@@ -38,6 +38,12 @@
namespace ceres {
namespace internal {
+// Noise factor for randomized cost function.
+static const double kNoiseFactor = 0.01;
+
+// Default random seed for randomized cost function.
+static const unsigned int kRandomSeed = 1234;
+
// y1 = x1'x2 -> dy1/dx1 = x2, dy1/dx2 = x1
// y2 = (x1'x2)^2 -> dy2/dx1 = 2 * x2 * (x1'x2), dy2/dx2 = 2 * x1 * (x1'x2)
// y3 = x2'x2 -> dy3/dx1 = 0, dy3/dx2 = 2 * x2
@@ -46,7 +52,7 @@
bool operator()(const double* x1, const double* x2, double* residuals) const;
void ExpectCostFunctionEvaluationIsNearlyCorrect(
const CostFunction& cost_function,
- NumericDiffMethod method) const;
+ NumericDiffMethodType method) const;
};
class EasyCostFunction : public SizedCostFunction<3, 5, 5> {
@@ -71,7 +77,7 @@
bool operator()(const double* x1, const double* x2, double* residuals) const;
void ExpectCostFunctionEvaluationIsNearlyCorrect(
const CostFunction& cost_function,
- NumericDiffMethod method) const;
+ NumericDiffMethodType method) const;
};
class TranscendentalCostFunction : public SizedCostFunction<2, 5, 5> {
@@ -85,6 +91,61 @@
TranscendentalFunctor functor_;
};
+// y = exp(x), dy/dx = exp(x)
+class ExponentialFunctor {
+ public:
+ bool operator()(const double* x1, double* residuals) const;
+ void ExpectCostFunctionEvaluationIsNearlyCorrect(
+ const CostFunction& cost_function) const;
+};
+
+class ExponentialCostFunction : public SizedCostFunction<1, 1> {
+ public:
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** /* not used */) const {
+ return functor_(parameters[0], residuals);
+ }
+
+ private:
+ ExponentialFunctor functor_;
+};
+
+// Test adaptive numeric differentiation by synthetically adding random noise
+// to a functor.
+// y = x^2 + [random noise], dy/dx ~ 2x
+class RandomizedFunctor {
+ public:
+ RandomizedFunctor(double noise_factor, unsigned int random_seed)
+ : noise_factor_(noise_factor), random_seed_(random_seed) {
+ }
+
+ bool operator()(const double* x1, double* residuals) const;
+ void ExpectCostFunctionEvaluationIsNearlyCorrect(
+ const CostFunction& cost_function) const;
+
+ private:
+ double noise_factor_;
+ unsigned int random_seed_;
+};
+
+class RandomizedCostFunction : public SizedCostFunction<1, 1> {
+ public:
+ RandomizedCostFunction(double noise_factor, unsigned int random_seed)
+ : functor_(noise_factor, random_seed) {
+ }
+
+ virtual bool Evaluate(double const* const* parameters,
+ double* residuals,
+ double** /* not used */) const {
+ return functor_(parameters[0], residuals);
+ }
+
+ private:
+ RandomizedFunctor functor_;
+};
+
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc
index 275d7f8..f86fb78 100644
--- a/internal/ceres/types.cc
+++ b/internal/ceres/types.cc
@@ -284,6 +284,27 @@
return false;
}
+const char* NumericDiffMethodTypeToString(
+ NumericDiffMethodType type) {
+ switch (type) {
+ CASESTR(CENTRAL);
+ CASESTR(FORWARD);
+ CASESTR(RIDDERS);
+ default:
+ return "UNKNOWN";
+ }
+}
+
+bool StringToNumericDiffMethodType(
+ string value,
+ NumericDiffMethodType* type) {
+ UpperCase(&value);
+ STRENUM(CENTRAL);
+ STRENUM(FORWARD);
+ STRENUM(RIDDERS);
+ return false;
+}
+
const char* VisibilityClusteringTypeToString(
VisibilityClusteringType type) {
switch (type) {
