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/docs/source/bibliography.rst b/docs/source/bibliography.rst
index bab6e61..ab2a29f 100644
--- a/docs/source/bibliography.rst
+++ b/docs/source/bibliography.rst
@@ -94,6 +94,9 @@
Part II: Implementation and Experiments**, Management Science,
20(5), 863-874, 1974.
+.. [Ridders] C. J. F. Ridders, **Accurate computation of F'(x) and
+ F'(x) F"(x)**, Advances in Engineering Software 4(2), 75-76, 1978.
+
.. [RuheWedin] A. Ruhe and P.Å. Wedin, **Algorithms for separable
nonlinear least squares problems**, Siam Review, 22(3):318–337,
1980.
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 408544d..038684d 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -367,7 +367,7 @@
.. code-block:: c++
template <typename CostFunctor,
- NumericDiffMethod method = CENTRAL,
+ NumericDiffMethodType method = CENTRAL,
int kNumResiduals, // Number of residuals, or ceres::DYNAMIC.
int N0, // Number of parameters in block 0.
int N1 = 0, // Number of parameters in block 1.
@@ -481,10 +481,30 @@
independent variables, and there is no limit on the dimensionality
of each of them.
- The ``CENTRAL`` difference method is considerably more accurate at
+ There are three available numeric differentiation schemes in ceres-solver:
+
+ The ``FORWARD`` difference method, which approximates :math:`f'(x)`
+ by computing :math:`\frac{f(x+h)-f(x)}{h}`, computes the cost function
+ one additional time at :math:`x+h`. It is the fastest but least accurate
+ method.
+
+ The ``CENTRAL`` difference method is more accurate at
the cost of twice as many function evaluations than forward
- difference. Consider using central differences begin with, and only
- after that works, trying forward difference to improve performance.
+ difference, estimating :math:`f'(x)` by computing
+ :math:`\frac{f(x+h)-f(x-h)}{2h}`.
+
+ The ``RIDDERS`` difference method[Ridders]_ is an adaptive scheme that
+ estimates derivatives by performing multiple central differences
+ at varying scales. Specifically, the algorithm starts at a certain
+ :math:`h` and as the derivative is estimated, this step size decreases.
+ To conserve function evaluations and estimate the derivative error, the
+ method performs Richardson extrapolations between the tested step sizes.
+ The algorithm exhibits considerably higher accuracy, but does so by
+ additional evaluations of the cost function.
+
+ Consider using ``CENTRAL`` differences to begin with. Based on the
+ results, either try forward difference to improve performance or
+ Ridders' method to improve accuracy.
**WARNING** A common beginner's error when first using
NumericDiffCostFunction is to get the sizing wrong. In particular,
@@ -530,7 +550,7 @@
**Alternate Interface**
- For a variety of reason, including compatibility with legacy code,
+ For a variety of reasons, including compatibility with legacy code,
:class:`NumericDiffCostFunction` can also take
:class:`CostFunction` objects as input. The following describes
how.
@@ -571,7 +591,7 @@
.. code-block:: c++
- template <typename CostFunctor, NumericDiffMethod method = CENTRAL>
+ template <typename CostFunctor, NumericDiffMethodType method = CENTRAL>
class DynamicNumericDiffCostFunction : public CostFunction {
};
diff --git a/examples/more_garbow_hillstrom.cc b/examples/more_garbow_hillstrom.cc
index b0ea78a..835cb7f 100644
--- a/examples/more_garbow_hillstrom.cc
+++ b/examples/more_garbow_hillstrom.cc
@@ -63,34 +63,62 @@
DEFINE_bool(use_numeric_diff, false,
"Use numeric differentiation instead of automatic "
"differentiation.");
-
+DEFINE_string(numeric_diff_method, "ridders", "When using numeric "
+ "differentiation, selects algorithm. Options are: central, "
+ "forward, ridders.");
+DEFINE_int32(ridders_extrapolations, 3, "Maximal number of extrapolations in "
+ "Ridders' method.");
namespace ceres {
namespace examples {
const double kDoubleMax = std::numeric_limits<double>::max();
-#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals) \
- struct name { \
- static const int kNumParameters = num_parameters; \
- static const double initial_x[kNumParameters]; \
- static const double lower_bounds[kNumParameters]; \
- static const double upper_bounds[kNumParameters]; \
- static const double constrained_optimal_cost; \
- static const double unconstrained_optimal_cost; \
- static CostFunction* Create() { \
- if (FLAGS_use_numeric_diff) { \
- return new NumericDiffCostFunction<name, \
- CENTRAL, \
- num_residuals, \
- num_parameters>(new name); \
- } else { \
- return new AutoDiffCostFunction<name, \
- num_residuals, \
- num_parameters>(new name); \
- } \
- } \
- template <typename T> \
+static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
+ options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
+}
+
+#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals) \
+ struct name { \
+ static const int kNumParameters = num_parameters; \
+ static const double initial_x[kNumParameters]; \
+ static const double lower_bounds[kNumParameters]; \
+ static const double upper_bounds[kNumParameters]; \
+ static const double constrained_optimal_cost; \
+ static const double unconstrained_optimal_cost; \
+ static CostFunction* Create() { \
+ if (FLAGS_use_numeric_diff) { \
+ ceres::NumericDiffOptions options; \
+ SetNumericDiffOptions(&options); \
+ if (FLAGS_numeric_diff_method == "central") { \
+ return new NumericDiffCostFunction<name, \
+ ceres::CENTRAL, \
+ num_residuals, \
+ num_parameters>( \
+ new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \
+ } else if (FLAGS_numeric_diff_method == "forward") { \
+ return new NumericDiffCostFunction<name, \
+ ceres::FORWARD, \
+ num_residuals, \
+ num_parameters>( \
+ new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \
+ } else if (FLAGS_numeric_diff_method == "ridders") { \
+ return new NumericDiffCostFunction<name, \
+ ceres::RIDDERS, \
+ num_residuals, \
+ num_parameters>( \
+ new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \
+ } else { \
+ LOG(ERROR) << "Invalid numeric diff method specified"; \
+ return NULL; \
+ } \
+ } else { \
+ return new AutoDiffCostFunction<name, \
+ num_residuals, \
+ num_parameters>(new name); \
+ } \
+ } \
+ template <typename T> \
bool operator()(const T* const x, T* residual) const {
#define END_MGH_PROBLEM return true; } }; // NOLINT
diff --git a/examples/nist.cc b/examples/nist.cc
index ceeb457..f7bd38e 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -65,7 +65,7 @@
// solver had the highest LRE.
// In this file, we implement the same evaluation methodology using
-// Ceres. Currently using Levenberg-Marquard with DENSE_QR, we get
+// Ceres. Currently using Levenberg-Marquardt with DENSE_QR, we get
//
// Excel Gnuplot GaussFit HBN MinPack Ceres
// Average LRE 2.3 4.3 4.0 6.8 4.4 9.4
@@ -118,6 +118,13 @@
DEFINE_bool(use_numeric_diff, false,
"Use numeric differentiation instead of automatic "
"differentiation.");
+DEFINE_string(numeric_diff_method, "ridders", "When using numeric "
+ "differentiation, selects algorithm. Options are: central, "
+ "forward, ridders.");
+DEFINE_double(ridders_step_size, 1e-9, "Initial step size for Ridders "
+ "numeric differentiation.");
+DEFINE_int32(ridders_extrapolations, 3, "Maximal number of Ridders "
+ "extrapolations.");
namespace ceres {
namespace examples {
@@ -413,6 +420,11 @@
double y_;
};
+static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
+ options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
+ options->ridders_relative_initial_step_size = FLAGS_ridders_step_size;
+}
+
template <typename Model, int num_residuals, int num_parameters>
int RegressionDriver(const string& filename,
const ceres::Solver::Options& options) {
@@ -439,11 +451,30 @@
response.data() + nist_problem.response_size() * i);
ceres::CostFunction* cost_function = NULL;
if (FLAGS_use_numeric_diff) {
- cost_function =
- new ceres::NumericDiffCostFunction<Model,
- ceres::CENTRAL,
- num_residuals,
- num_parameters>(model);
+ ceres::NumericDiffOptions options;
+ SetNumericDiffOptions(&options);
+ if (FLAGS_numeric_diff_method == "central") {
+ cost_function = new NumericDiffCostFunction<Model,
+ ceres::CENTRAL,
+ num_residuals,
+ num_parameters>(
+ model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+ } else if (FLAGS_numeric_diff_method == "forward") {
+ cost_function = new NumericDiffCostFunction<Model,
+ ceres::FORWARD,
+ num_residuals,
+ num_parameters>(
+ model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+ } else if (FLAGS_numeric_diff_method == "ridders") {
+ cost_function = new NumericDiffCostFunction<Model,
+ ceres::RIDDERS,
+ num_residuals,
+ num_parameters>(
+ model, ceres::TAKE_OWNERSHIP, num_residuals, options);
+ } else {
+ LOG(ERROR) << "Invalid numeric diff method specified";
+ return 0;
+ }
} else {
cost_function =
new ceres::AutoDiffCostFunction<Model,
diff --git a/include/ceres/dynamic_numeric_diff_cost_function.h b/include/ceres/dynamic_numeric_diff_cost_function.h
index 72fa13b..c852d57 100644
--- a/include/ceres/dynamic_numeric_diff_cost_function.h
+++ b/include/ceres/dynamic_numeric_diff_cost_function.h
@@ -68,19 +68,37 @@
#include "ceres/internal/scoped_ptr.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/numeric_diff.h"
+#include "ceres/numeric_diff_options.h"
#include "glog/logging.h"
namespace ceres {
-template <typename CostFunctor, NumericDiffMethod method = CENTRAL>
+template <typename CostFunctor, NumericDiffMethodType method = CENTRAL>
class DynamicNumericDiffCostFunction : public CostFunction {
public:
- explicit DynamicNumericDiffCostFunction(const CostFunctor* functor,
- Ownership ownership = TAKE_OWNERSHIP,
- double relative_step_size = 1e-6)
+ explicit DynamicNumericDiffCostFunction(
+ const CostFunctor* functor,
+ Ownership ownership = TAKE_OWNERSHIP,
+ const NumericDiffOptions& options = NumericDiffOptions())
: functor_(functor),
ownership_(ownership),
- relative_step_size_(relative_step_size) {
+ options_(options) {
+ }
+
+ // Deprecated. New users should avoid using this constructor. Instead, use the
+ // constructor with NumericDiffOptions.
+ DynamicNumericDiffCostFunction(
+ const CostFunctor* functor,
+ Ownership ownership,
+ double relative_step_size)
+ : functor_(functor),
+ ownership_(ownership),
+ options_() {
+ LOG(WARNING) << "This constructor is deprecated and will be removed in "
+ "a future version. Please use the NumericDiffOptions "
+ "constructor instead.";
+
+ options_.relative_step_size = relative_step_size;
}
virtual ~DynamicNumericDiffCostFunction() {
@@ -140,7 +158,7 @@
DYNAMIC, DYNAMIC>::EvaluateJacobianForParameterBlock(
functor_.get(),
residuals,
- relative_step_size_,
+ options_,
this->num_residuals(),
block,
block_sizes[block],
@@ -179,7 +197,7 @@
internal::scoped_ptr<const CostFunctor> functor_;
Ownership ownership_;
- const double relative_step_size_;
+ NumericDiffOptions options_;
};
} // namespace ceres
diff --git a/include/ceres/internal/numeric_diff.h b/include/ceres/internal/numeric_diff.h
index b9aab88..2fffdff 100644
--- a/include/ceres/internal/numeric_diff.h
+++ b/include/ceres/internal/numeric_diff.h
@@ -28,8 +28,9 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
// mierle@gmail.com (Keir Mierle)
+// tbennun@gmail.com (Tal Ben-Nun)
//
-// Finite differencing routine used by NumericDiffCostFunction.
+// Finite differencing routines used by NumericDiffCostFunction.
#ifndef CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
#define CERES_PUBLIC_INTERNAL_NUMERIC_DIFF_H_
@@ -37,9 +38,12 @@
#include <cstring>
#include "Eigen/Dense"
+#include "Eigen/StdVector"
#include "ceres/cost_function.h"
+#include "ceres/internal/fixed_array.h"
#include "ceres/internal/scoped_ptr.h"
#include "ceres/internal/variadic_evaluate.h"
+#include "ceres/numeric_diff_options.h"
#include "ceres/types.h"
#include "glog/logging.h"
@@ -78,7 +82,7 @@
// specializations for member functions. The alternative is to repeat the main
// class for differing numbers of parameters, which is also unfortunate.
template <typename CostFunctor,
- NumericDiffMethod kMethod,
+ NumericDiffMethodType kMethod,
int kNumResiduals,
int N0, int N1, int N2, int N3, int N4,
int N5, int N6, int N7, int N8, int N9,
@@ -88,8 +92,8 @@
// Mutates parameters but must restore them before return.
static bool EvaluateJacobianForParameterBlock(
const CostFunctor* functor,
- double const* residuals_at_eval_point,
- const double relative_step_size,
+ const double* residuals_at_eval_point,
+ const NumericDiffOptions& options,
int num_residuals,
int parameter_block_index,
int parameter_block_size,
@@ -126,67 +130,285 @@
num_residuals_internal,
parameter_block_size_internal);
- // Mutate 1 element at a time and then restore.
Map<ParameterVector> x_plus_delta(
parameters[parameter_block_index_internal],
parameter_block_size_internal);
ParameterVector x(x_plus_delta);
- ParameterVector step_size = x.array().abs() * relative_step_size;
+ ParameterVector step_size = x.array().abs() *
+ ((kMethod == RIDDERS) ? options.ridders_relative_initial_step_size :
+ options.relative_step_size);
// 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());
+ double min_step_size = std::sqrt(std::numeric_limits<double>::epsilon());
+
+ // For Ridders' method, the initial step size is required to be large,
+ // thus ridders_relative_initial_step_size is used.
+ if (kMethod == RIDDERS) {
+ min_step_size = std::max(min_step_size,
+ options.ridders_relative_initial_step_size);
+ }
// For each parameter in the parameter block, use finite differences to
// compute the derivative for that parameter.
- ResidualVector residuals(num_residuals_internal);
+ FixedArray<double> temp_residual_array(num_residuals_internal);
+ FixedArray<double> residual_array(num_residuals_internal);
+ Map<ResidualVector> residuals(residual_array.get(),
+ num_residuals_internal);
+
for (int j = 0; j < parameter_block_size_internal; ++j) {
const double delta = std::max(min_step_size, step_size(j));
- x_plus_delta(j) = x(j) + delta;
+
+ if (kMethod == RIDDERS) {
+ if (!EvaluateRiddersJacobianColumn(functor, j, delta,
+ options,
+ num_residuals_internal,
+ parameter_block_size_internal,
+ x.data(),
+ residuals_at_eval_point,
+ parameters,
+ x_plus_delta.data(),
+ temp_residual_array.get(),
+ residual_array.get())) {
+ return false;
+ }
+ } else {
+ if (!EvaluateJacobianColumn(functor, j, delta,
+ num_residuals_internal,
+ parameter_block_size_internal,
+ x.data(),
+ residuals_at_eval_point,
+ parameters,
+ x_plus_delta.data(),
+ temp_residual_array.get(),
+ residual_array.get())) {
+ return false;
+ }
+ }
+
+ parameter_jacobian.col(j).matrix() = residuals;
+ }
+ return true;
+ }
+
+ static bool EvaluateJacobianColumn(const CostFunctor* functor,
+ int parameter_index,
+ double delta,
+ int num_residuals,
+ int parameter_block_size,
+ const double* x_ptr,
+ const double* residuals_at_eval_point,
+ double** parameters,
+ double* x_plus_delta_ptr,
+ double* temp_residuals_ptr,
+ double* residuals_ptr) {
+ using Eigen::Map;
+ using Eigen::Matrix;
+
+ typedef Matrix<double, kNumResiduals, 1> ResidualVector;
+ typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
+
+ Map<const ParameterVector> x(x_ptr, parameter_block_size);
+ Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
+ parameter_block_size);
+
+ Map<ResidualVector> residuals(residuals_ptr, num_residuals);
+ Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
+
+ // Mutate 1 element at a time and then restore.
+ x_plus_delta(parameter_index) = x(parameter_index) + delta;
+
+ if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
+ functor, parameters, residuals.data(), functor)) {
+ return false;
+ }
+
+ // Compute this column of the jacobian in 3 steps:
+ // 1. Store residuals for the forward part.
+ // 2. Subtract residuals for the backward (or 0) part.
+ // 3. Divide out the run.
+ double one_over_delta = 1.0 / delta;
+ if (kMethod == CENTRAL || kMethod == RIDDERS) {
+ // Compute the function on the other side of x(parameter_index).
+ x_plus_delta(parameter_index) = x(parameter_index) - delta;
if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
- functor, parameters, residuals.data(), functor)) {
+ functor, parameters, temp_residuals.data(), functor)) {
return false;
}
- // Compute this column of the jacobian in 3 steps:
- // 1. Store residuals for the forward part.
- // 2. Subtract residuals for the backward (or 0) part.
- // 3. Divide out the run.
- parameter_jacobian.col(j) = residuals;
+ residuals -= temp_residuals;
+ one_over_delta /= 2;
+ } else {
+ // Forward difference only; reuse existing residuals evaluation.
+ residuals -=
+ Map<const ResidualVector>(residuals_at_eval_point,
+ num_residuals);
+ }
- double one_over_delta = 1.0 / delta;
- if (kMethod == CENTRAL) {
- // Compute the function on the other side of x(j).
- x_plus_delta(j) = x(j) - delta;
+ // Restore x_plus_delta.
+ x_plus_delta(parameter_index) = x(parameter_index);
- if (!EvaluateImpl<CostFunctor, N0, N1, N2, N3, N4, N5, N6, N7, N8, N9>(
- functor, parameters, residuals.data(), functor)) {
- return false;
- }
+ // Divide out the run to get slope.
+ residuals *= one_over_delta;
- parameter_jacobian.col(j) -= residuals;
- one_over_delta /= 2;
- } else {
- // Forward difference only; reuse existing residuals evaluation.
- parameter_jacobian.col(j) -=
- Map<const ResidualVector>(residuals_at_eval_point,
- num_residuals_internal);
+ return true;
+ }
+
+ // This numeric difference implementation uses adaptive differentiation
+ // on the parameters to obtain the Jacobian matrix. The adaptive algorithm
+ // is based on Ridders' method for adaptive differentiation, which creates
+ // a Romberg tableau from varying step sizes and extrapolates the
+ // intermediate results to obtain the current computational error.
+ //
+ // References:
+ // C.J.F. Ridders, Accurate computation of F'(x) and F'(x) F"(x), Advances
+ // in Engineering Software (1978), Volume 4, Issue 2, April 1982,
+ // Pages 75-76, ISSN 0141-1195,
+ // http://dx.doi.org/10.1016/S0141-1195(82)80057-0.
+ static bool EvaluateRiddersJacobianColumn(
+ const CostFunctor* functor,
+ int parameter_index,
+ double delta,
+ const NumericDiffOptions& options,
+ int num_residuals,
+ int parameter_block_size,
+ const double* x_ptr,
+ const double* residuals_at_eval_point,
+ double** parameters,
+ double* x_plus_delta_ptr,
+ double* temp_residuals_ptr,
+ double* residuals_ptr) {
+ using Eigen::Map;
+ using Eigen::Matrix;
+ using Eigen::aligned_allocator;
+
+ typedef Matrix<double, kNumResiduals, 1> ResidualVector;
+ typedef Matrix<double, kNumResiduals, DYNAMIC> ResidualCandidateMatrix;
+ typedef Matrix<double, kParameterBlockSize, 1> ParameterVector;
+
+ Map<const ParameterVector> x(x_ptr, parameter_block_size);
+ Map<ParameterVector> x_plus_delta(x_plus_delta_ptr,
+ parameter_block_size);
+
+ Map<ResidualVector> residuals(residuals_ptr, num_residuals);
+ Map<ResidualVector> temp_residuals(temp_residuals_ptr, num_residuals);
+
+ // In order for the algorithm to converge, the step size should be
+ // initialized to a value that is large enough to produce a significant
+ // change in the function.
+ // As the derivative is estimated, the step size decreases.
+ // By default, the step sizes are chosen so that the middle column
+ // of the Romberg tableau uses the input delta.
+ double current_step_size = delta *
+ pow(options.ridders_step_shrink_factor,
+ options.max_num_ridders_extrapolations / 2);
+
+ // Double-buffering temporary differential candidate vectors
+ // from previous step size.
+ ResidualCandidateMatrix stepsize_candidates_a(
+ num_residuals,
+ options.max_num_ridders_extrapolations);
+ ResidualCandidateMatrix stepsize_candidates_b(
+ num_residuals,
+ options.max_num_ridders_extrapolations);
+ ResidualCandidateMatrix* current_candidates = &stepsize_candidates_a;
+ ResidualCandidateMatrix* previous_candidates = &stepsize_candidates_b;
+
+ // Represents the computational error of the derivative. This variable is
+ // initially set to a large value, and is set to the difference between
+ // current and previous finite difference extrapolations.
+ // norm_error is supposed to decrease as the finite difference tableau
+ // generation progresses, serving both as an estimate for differentiation
+ // error and as a measure of differentiation numerical stability.
+ double norm_error = std::numeric_limits<double>::max();
+
+ // Loop over decreasing step sizes until:
+ // 1. Error is smaller than a given value (ridders_epsilon),
+ // 2. Maximal order of extrapolation reached, or
+ // 3. Extrapolation becomes numerically unstable.
+ for (int i = 0; i < options.max_num_ridders_extrapolations; ++i) {
+ // Compute the numerical derivative at this step size.
+ if (!EvaluateJacobianColumn(functor, parameter_index, current_step_size,
+ num_residuals,
+ parameter_block_size,
+ x.data(),
+ residuals_at_eval_point,
+ parameters,
+ x_plus_delta.data(),
+ temp_residuals.data(),
+ current_candidates->col(0).data())) {
+ // Something went wrong; bail.
+ return false;
}
- x_plus_delta(j) = x(j); // Restore x_plus_delta.
- // Divide out the run to get slope.
- parameter_jacobian.col(j) *= one_over_delta;
+ // Store initial results.
+ if (i == 0) {
+ residuals = current_candidates[0];
+ }
+
+ // Shrink differentiation step size.
+ current_step_size /= options.ridders_step_shrink_factor;
+
+ // Extrapolation factor for Richardson acceleration method (see below).
+ double richardson_factor = options.ridders_step_shrink_factor *
+ options.ridders_step_shrink_factor;
+ for (int k = 1; k <= i; ++k) {
+ // Extrapolate the various orders of finite differences using
+ // the Richardson acceleration method.
+ current_candidates->col(k) =
+ (richardson_factor * current_candidates->col(k - 1) -
+ previous_candidates->col(k - 1)) / (richardson_factor - 1.0);
+
+ richardson_factor *= options.ridders_step_shrink_factor *
+ options.ridders_step_shrink_factor;
+
+ // Compute the difference between the previous value and the current.
+ double candidate_error = std::max(
+ (current_candidates->col(k) -
+ current_candidates->col(k - 1)).norm(),
+ (current_candidates->col(k) -
+ previous_candidates->col(k - 1)).norm());
+
+ // If the error has decreased, update results.
+ if (candidate_error <= norm_error) {
+ norm_error = candidate_error;
+ residuals = current_candidates->col(k);
+
+ // If the error is small enough, stop.
+ if (norm_error < options.ridders_epsilon) {
+ break;
+ }
+ }
+ }
+
+ // After breaking out of the inner loop, declare convergence.
+ if (norm_error < options.ridders_epsilon) {
+ break;
+ }
+
+ // Check to see if the current gradient estimate is numerically unstable.
+ // If so, bail out and return the last stable result.
+ if (i > 0) {
+ double tableau_error = (current_candidates->col(i) -
+ previous_candidates->col(i - 1)).norm();
+
+ // Compare current error to the chosen candidate's error.
+ if (tableau_error >= 2 * norm_error) {
+ break;
+ }
+ }
+
+ std::swap(current_candidates, previous_candidates);
}
return true;
}
};
template <typename CostFunctor,
- NumericDiffMethod kMethod,
+ NumericDiffMethodType kMethod,
int kNumResiduals,
int N0, int N1, int N2, int N3, int N4,
int N5, int N6, int N7, int N8, int N9,
@@ -197,8 +419,8 @@
// Mutates parameters but must restore them before return.
static bool EvaluateJacobianForParameterBlock(
const CostFunctor* functor,
- double const* residuals_at_eval_point,
- const double relative_step_size,
+ const double* residuals_at_eval_point,
+ const NumericDiffOptions& options,
const int num_residuals,
const int parameter_block_index,
const int parameter_block_size,
@@ -207,7 +429,7 @@
// Silence unused parameter compiler warnings.
(void)functor;
(void)residuals_at_eval_point;
- (void)relative_step_size;
+ (void)options;
(void)num_residuals;
(void)parameter_block_index;
(void)parameter_block_size;
diff --git a/include/ceres/numeric_diff_cost_function.h b/include/ceres/numeric_diff_cost_function.h
index c15e06d..fa96078 100644
--- a/include/ceres/numeric_diff_cost_function.h
+++ b/include/ceres/numeric_diff_cost_function.h
@@ -132,7 +132,7 @@
//
// ALTERNATE INTERFACE
//
-// For a variety of reason, including compatibility with legacy code,
+// For a variety of reasons, including compatibility with legacy code,
// NumericDiffCostFunction can also take CostFunction objects as
// input. The following describes how.
//
@@ -165,6 +165,7 @@
#include "ceres/cost_function.h"
#include "ceres/internal/numeric_diff.h"
#include "ceres/internal/scoped_ptr.h"
+#include "ceres/numeric_diff_options.h"
#include "ceres/sized_cost_function.h"
#include "ceres/types.h"
#include "glog/logging.h"
@@ -172,7 +173,7 @@
namespace ceres {
template <typename CostFunctor,
- NumericDiffMethod method = CENTRAL,
+ NumericDiffMethodType method = CENTRAL,
int kNumResiduals = 0, // Number of residuals, or ceres::DYNAMIC
int N0 = 0, // Number of parameters in block 0.
int N1 = 0, // Number of parameters in block 1.
@@ -189,13 +190,14 @@
N0, N1, N2, N3, N4,
N5, N6, N7, N8, N9> {
public:
- NumericDiffCostFunction(CostFunctor* functor,
- Ownership ownership = TAKE_OWNERSHIP,
- int num_residuals = kNumResiduals,
- const double relative_step_size = 1e-6)
- :functor_(functor),
- ownership_(ownership),
- relative_step_size_(relative_step_size) {
+ NumericDiffCostFunction(
+ CostFunctor* functor,
+ Ownership ownership = TAKE_OWNERSHIP,
+ int num_residuals = kNumResiduals,
+ const NumericDiffOptions& options = NumericDiffOptions())
+ : functor_(functor),
+ ownership_(ownership),
+ options_(options) {
if (kNumResiduals == DYNAMIC) {
SizedCostFunction<kNumResiduals,
N0, N1, N2, N3, N4,
@@ -204,6 +206,29 @@
}
}
+ // Deprecated. New users should avoid using this constructor. Instead, use the
+ // constructor with NumericDiffOptions.
+ NumericDiffCostFunction(CostFunctor* functor,
+ Ownership ownership,
+ int num_residuals,
+ const double relative_step_size)
+ :functor_(functor),
+ ownership_(ownership),
+ options_() {
+ LOG(WARNING) << "This constructor is deprecated and will be removed in "
+ "a future version. Please use the NumericDiffOptions "
+ "constructor instead.";
+
+ if (kNumResiduals == DYNAMIC) {
+ SizedCostFunction<kNumResiduals,
+ N0, N1, N2, N3, N4,
+ N5, N6, N7, N8, N9>
+ ::set_num_residuals(num_residuals);
+ }
+
+ options_.relative_step_size = relative_step_size;
+ }
+
~NumericDiffCostFunction() {
if (ownership_ != TAKE_OWNERSHIP) {
functor_.release();
@@ -278,7 +303,7 @@
N ## block >::EvaluateJacobianForParameterBlock( \
functor_.get(), \
residuals, \
- relative_step_size_, \
+ options_, \
SizedCostFunction<kNumResiduals, \
N0, N1, N2, N3, N4, \
N5, N6, N7, N8, N9>::num_residuals(), \
@@ -309,7 +334,7 @@
private:
internal::scoped_ptr<CostFunctor> functor_;
Ownership ownership_;
- const double relative_step_size_;
+ NumericDiffOptions options_;
};
} // namespace ceres
diff --git a/include/ceres/numeric_diff_options.h b/include/ceres/numeric_diff_options.h
new file mode 100644
index 0000000..119c8a8
--- /dev/null
+++ b/include/ceres/numeric_diff_options.h
@@ -0,0 +1,79 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2015 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: tbennun@gmail.com (Tal Ben-Nun)
+//
+
+#ifndef CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
+#define CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
+
+namespace ceres {
+
+// Options pertaining to numeric differentiation (e.g., convergence criteria,
+// step sizes).
+struct CERES_EXPORT NumericDiffOptions {
+ NumericDiffOptions() {
+ relative_step_size = 1e-6;
+ ridders_relative_initial_step_size = 1e-2;
+ max_num_ridders_extrapolations = 10;
+ ridders_epsilon = 1e-12;
+ ridders_step_shrink_factor = 2.0;
+ }
+
+ // Numeric differentiation step size (multiplied by parameter block's
+ // order of magnitude). If parameters are close to zero, the step size
+ // is set to sqrt(machine_epsilon).
+ double relative_step_size;
+
+ // Initial step size for Ridders adaptive numeric differentiation (multiplied
+ // by parameter block's order of magnitude).
+ // If parameters are close to zero, Ridders' method sets the step size
+ // directly to this value. This parameter is separate from
+ // "relative_step_size" in order to set a different default value.
+ //
+ // Note: For Ridders' method to converge, the step size should be initialized
+ // to a value that is large enough to produce a significant change in the
+ // function. As the derivative is estimated, the step size decreases.
+ double ridders_relative_initial_step_size;
+
+ // Maximal number of adaptive extrapolations (sampling) in Ridders' method.
+ int max_num_ridders_extrapolations;
+
+ // Convergence criterion on extrapolation error for Ridders adaptive
+ // differentiation. The available error estimation methods are defined in
+ // NumericDiffErrorType and set in the "ridders_error_method" field.
+ double ridders_epsilon;
+
+ // The factor in which to shrink the step size with each extrapolation in
+ // Ridders' method.
+ double ridders_step_shrink_factor;
+};
+
+} // namespace ceres
+
+#endif // CERES_PUBLIC_NUMERIC_DIFF_OPTIONS_H_
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 9b01af2..2ea4180 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -397,9 +397,19 @@
DYNAMIC = -1
};
-enum NumericDiffMethod {
+// The differentiation method used to compute numerical derivatives in
+// NumericDiffCostFunction and DynamicNumericDiffCostFunction.
+enum NumericDiffMethodType {
+ // Compute central finite difference: f'(x) ~ (f(x+h) - f(x-h)) / 2h.
CENTRAL,
- FORWARD
+
+ // Compute forward finite difference: f'(x) ~ (f(x+h) - f(x)) / h.
+ FORWARD,
+
+ // Adaptive numerical differentiation using Ridders' method. Provides more
+ // accurate and robust derivatives at the expense of additional cost
+ // function evaluations.
+ RIDDERS
};
enum LineSearchInterpolationType {
@@ -477,6 +487,12 @@
std::string value,
CovarianceAlgorithmType* type);
+CERES_EXPORT const char* NumericDiffMethodTypeToString(
+ NumericDiffMethodType type);
+CERES_EXPORT bool StringToNumericDiffMethodType(
+ std::string value,
+ NumericDiffMethodType* type);
+
CERES_EXPORT const char* TerminationTypeToString(TerminationType type);
CERES_EXPORT bool IsSchurType(LinearSolverType type);
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) {
