diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst index bab6e6171..ab2a29f71 100644 --- a/docs/source/bibliography.rst +++ b/docs/source/bibliography.rst @@ -94,6 +94,9 @@ Bibliography 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 408544d51..038684d10 100644 --- a/docs/source/nnls_modeling.rst +++ b/docs/source/nnls_modeling.rst @@ -367,7 +367,7 @@ the corresponding accessors. This information will be verified by the .. code-block:: c++ template + template class DynamicNumericDiffCostFunction : public CostFunction { }; diff --git a/examples/more_garbow_hillstrom.cc b/examples/more_garbow_hillstrom.cc index b0ea78a99..835cb7fd6 100644 --- a/examples/more_garbow_hillstrom.cc +++ b/examples/more_garbow_hillstrom.cc @@ -63,34 +63,62 @@ DEFINE_string(problem, "all", "Which problem to solve"); 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::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(new name); \ - } else { \ - return new AutoDiffCostFunction(new name); \ - } \ - } \ - template \ +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( \ + new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \ + } else if (FLAGS_numeric_diff_method == "forward") { \ + return new NumericDiffCostFunction( \ + new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \ + } else if (FLAGS_numeric_diff_method == "ridders") { \ + return new NumericDiffCostFunction( \ + new name, ceres::TAKE_OWNERSHIP, num_residuals, options); \ + } else { \ + LOG(ERROR) << "Invalid numeric diff method specified"; \ + return NULL; \ + } \ + } else { \ + return new AutoDiffCostFunction(new name); \ + } \ + } \ + template \ 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 ceeb4578d..f7bd38e1f 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_double(initial_trust_region_radius, 1e4, "Initial trust region radius"); 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 @@ struct Nelson { 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 int RegressionDriver(const string& filename, const ceres::Solver::Options& options) { @@ -439,11 +451,30 @@ int RegressionDriver(const string& filename, response.data() + nist_problem.response_size() * i); ceres::CostFunction* cost_function = NULL; if (FLAGS_use_numeric_diff) { - cost_function = - new ceres::NumericDiffCostFunction(model); + ceres::NumericDiffOptions options; + SetNumericDiffOptions(&options); + if (FLAGS_numeric_diff_method == "central") { + cost_function = new NumericDiffCostFunction( + model, ceres::TAKE_OWNERSHIP, num_residuals, options); + } else if (FLAGS_numeric_diff_method == "forward") { + cost_function = new NumericDiffCostFunction( + model, ceres::TAKE_OWNERSHIP, num_residuals, options); + } else if (FLAGS_numeric_diff_method == "ridders") { + cost_function = new NumericDiffCostFunction( + model, ceres::TAKE_OWNERSHIP, num_residuals, options); + } else { + LOG(ERROR) << "Invalid numeric diff method specified"; + return 0; + } } else { cost_function = new ceres::AutoDiffCostFunction +template 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 @@ class DynamicNumericDiffCostFunction : public CostFunction { DYNAMIC, DYNAMIC>::EvaluateJacobianForParameterBlock( functor_.get(), residuals, - relative_step_size_, + options_, this->num_residuals(), block, block_sizes[block], @@ -179,7 +197,7 @@ class DynamicNumericDiffCostFunction : public CostFunction { internal::scoped_ptr 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 b9aab88e6..2fffdffe6 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 #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 @@ bool EvaluateImpl(const CostFunctor* functor, // specializations for member functions. The alternative is to repeat the main // class for differing numbers of parameters, which is also unfortunate. template 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::epsilon()); + double min_step_size = std::sqrt(std::numeric_limits::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 temp_residual_array(num_residuals_internal); + FixedArray residual_array(num_residuals_internal); + Map 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 ResidualVector; + typedef Matrix ParameterVector; + + Map x(x_ptr, parameter_block_size); + Map x_plus_delta(x_plus_delta_ptr, + parameter_block_size); + + Map residuals(residuals_ptr, num_residuals); + Map 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( + 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( - 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(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( - 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(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 ResidualVector; + typedef Matrix ResidualCandidateMatrix; + typedef Matrix ParameterVector; + + Map x(x_ptr, parameter_block_size); + Map x_plus_delta(x_plus_delta_ptr, + parameter_block_size); + + Map residuals(residuals_ptr, num_residuals); + Map 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::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 { 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 + ::set_num_residuals(num_residuals); + } + + options_.relative_step_size = relative_step_size; + } + ~NumericDiffCostFunction() { if (ownership_ != TAKE_OWNERSHIP) { functor_.release(); @@ -278,7 +303,7 @@ class NumericDiffCostFunction N ## block >::EvaluateJacobianForParameterBlock( \ functor_.get(), \ residuals, \ - relative_step_size_, \ + options_, \ SizedCostFunction::num_residuals(), \ @@ -309,7 +334,7 @@ class NumericDiffCostFunction private: internal::scoped_ptr 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 000000000..119c8a865 --- /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 9b01af2d5..2ea418036 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -397,9 +397,19 @@ enum DimensionType { 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 @@ CERES_EXPORT bool StringToCovarianceAlgorithmType( 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 e3ad3d459..580fd260e 100644 --- a/internal/ceres/gradient_checking_cost_function.cc +++ b/internal/ceres/gradient_checking_cost_function.cc @@ -86,7 +86,7 @@ bool IsClose(double x, double y, double relative_precision, 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 @@ class GradientCheckingCostFunction : public CostFunction { new DynamicNumericDiffCostFunction( function, DO_NOT_TAKE_OWNERSHIP, - relative_step_size); + options); const vector& parameter_block_sizes = function->parameter_block_sizes(); @@ -235,8 +235,11 @@ CostFunction *CreateGradientCheckingCostFunction( 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 cb74a86a3..13ab1069c 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 @@ TEST(NumericDiffCostFunction, EasyCaseFunctorForwardDifferences) { functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD); } +TEST(NumericDiffCostFunction, EasyCaseFunctorRidders) { + internal::scoped_ptr cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new EasyFunctor)); + EasyFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS); +} + TEST(NumericDiffCostFunction, EasyCaseCostFunctionCentralDifferences) { internal::scoped_ptr cost_function; cost_function.reset( @@ -97,7 +111,21 @@ TEST(NumericDiffCostFunction, EasyCaseCostFunctionForwardDifferences) { functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, FORWARD); } -TEST(NumericDiffCostFunction, TranscendentalCaseFunctorCentralDifferences) { +TEST(NumericDiffCostFunction, EasyCaseCostFunctionRidders) { + internal::scoped_ptr cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new EasyCostFunction, TAKE_OWNERSHIP)); + EasyFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS); +} + +TEST(NumericDiffCostFunction, + TranscendentalCaseFunctorCentralDifferences) { internal::scoped_ptr cost_function; cost_function.reset( new NumericDiffCostFunction cost_function; cost_function.reset( new NumericDiffCostFunction cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new TranscendentalFunctor, TAKE_OWNERSHIP, 2, options)); + TranscendentalFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS); +} + +TEST(NumericDiffCostFunction, + TranscendentalCaseCostFunctionCentralDifferences) { internal::scoped_ptr cost_function; cost_function.reset( new NumericDiffCostFunction cost_function; cost_function.reset( new NumericDiffCostFunction cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new TranscendentalCostFunction, TAKE_OWNERSHIP, 2, options)); + TranscendentalFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function, RIDDERS); +} + template class SizeTestingCostFunction : public SizedCostFunction { public: @@ -208,7 +279,8 @@ TEST(NumericDiffCostFunction, EigenRowMajorColMajorTest) { new EasyFunctor, TAKE_OWNERSHIP, 2)); } -TEST(NumericDiffCostFunction, EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) { +TEST(NumericDiffCostFunction, + EasyCaseFunctorCentralDifferencesAndDynamicNumResiduals) { internal::scoped_ptr cost_function; cost_function.reset( new NumericDiffCostFunction cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new ExponentialFunctor)); + ExponentialFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function); +} + +TEST(NumericDiffCostFunction, ExponentialCostFunctionRidders) { + internal::scoped_ptr cost_function; + cost_function.reset( + new NumericDiffCostFunction( + new ExponentialCostFunction)); + ExponentialFunctor functor; + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function); +} + +TEST(NumericDiffCostFunction, RandomizedFunctorRidders) { + internal::scoped_ptr 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( + new RandomizedFunctor(kNoiseFactor, kRandomSeed), TAKE_OWNERSHIP, + 1, options)); + RandomizedFunctor functor (kNoiseFactor, kRandomSeed); + functor.ExpectCostFunctionEvaluationIsNearlyCorrect(*cost_function); +} + +TEST(NumericDiffCostFunction, RandomizedCostFunctionRidders) { + internal::scoped_ptr 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( + 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 55954497c..bab186cb1 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 @@ bool EasyFunctor::operator()(const double* x1, 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 @@ void EasyFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect( 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 @@ bool TranscendentalFunctor::operator()(const double* x1, void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect( const CostFunction& cost_function, - NumericDiffMethod method) const { + NumericDiffMethodType method) const { struct { double x1[5]; double x2[5]; @@ -150,7 +165,21 @@ void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect( 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 @@ void TranscendentalFunctor::ExpectCostFunctionEvaluationIsNearlyCorrect( } } +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(rand()) / + static_cast(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 3ae83faf0..2a551d394 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 @@ class EasyFunctor { 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 @@ class TranscendentalFunctor { 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 @@ class TranscendentalCostFunction : public SizedCostFunction<2, 5, 5> { 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 275d7f8b2..f86fb78eb 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -284,6 +284,27 @@ bool StringToCovarianceAlgorithmType( 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) {