Refactoring of the LineSearchMinimizer. 1. New LineSearchDirection interface, factory and instances. 2. Cleanup of LineSearchMinimizer to use the State and Direction objects. 3. LBFGS -> LowRankInverseHessian. 4. Refactoring of the RunCallbacks function and share it across LineSearchMinimizer and TrustRegionMinimizer. Change-Id: I19354afc6f5d6567b28918710c2012dc30ef8f32
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 6b664e4..786a3d7 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -57,15 +57,17 @@ gradient_checking_cost_function.cc implicit_schur_complement.cc iterative_schur_complement_solver.cc - lbfgs.cc levenberg_marquardt_strategy.cc line_search.cc + line_search_direction.cc line_search_minimizer.cc linear_least_squares_problems.cc linear_operator.cc linear_solver.cc local_parameterization.cc loss_function.cc + low_rank_inverse_hessian.cc + minimizer.cc normal_prior.cc parameter_block_ordering.cc partitioned_matrix_view.cc
diff --git a/internal/ceres/line_search_direction.cc b/internal/ceres/line_search_direction.cc new file mode 100644 index 0000000..f441607 --- /dev/null +++ b/internal/ceres/line_search_direction.cc
@@ -0,0 +1,144 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2012 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/line_search_direction.h" +#include "ceres/line_search_minimizer.h" +#include "ceres/low_rank_inverse_hessian.h" +#include "ceres/internal/eigen.h" +#include "glog/logging.h" + +namespace ceres { +namespace internal { + +class SteepestDescent : public LineSearchDirection { + public: + virtual ~SteepestDescent() {} + bool NextDirection(const LineSearchMinimizer::State& previous, + const LineSearchMinimizer::State& current, + Vector* search_direction) { + *search_direction = -current.gradient; + return true; + } +}; + +class NonlinearConjugateGradient : public LineSearchDirection { + public: + NonlinearConjugateGradient(const NonlinearConjugateGradientType type, + const double function_tolerance) + : type_(type), + function_tolerance_(function_tolerance) { + } + + bool NextDirection(const LineSearchMinimizer::State& previous, + const LineSearchMinimizer::State& current, + Vector* search_direction) { + double beta = 0.0; + Vector gradient_change; + switch (type_) { + case FLETCHER_REEVES: + beta = current.gradient_squared_norm / previous.gradient_squared_norm; + break; + case POLAK_RIBIRERE: + gradient_change = current.gradient - previous.gradient; + beta = (current.gradient.dot(gradient_change) / + previous.gradient_squared_norm); + break; + case HESTENES_STIEFEL: + gradient_change = current.gradient - previous.gradient; + beta = (current.gradient.dot(gradient_change) / + previous.search_direction.dot(gradient_change)); + break; + default: + LOG(FATAL) << "Unknown nonlinear conjugate gradient type: " << type_; + } + + *search_direction = -current.gradient + beta * previous.search_direction; + const double directional_derivative = current. gradient.dot(*search_direction); + if (directional_derivative > -function_tolerance_) { + LOG(WARNING) << "Restarting non-linear conjugate gradients: " + << directional_derivative; + *search_direction = -current.gradient; + }; + + return true; + } + + private: + const NonlinearConjugateGradientType type_; + const double function_tolerance_; +}; + +class LBFGS : public LineSearchDirection { + public: + LBFGS(const int num_parameters, const int max_lbfgs_rank) + : low_rank_inverse_hessian_(num_parameters, max_lbfgs_rank) {} + + virtual ~LBFGS() {} + + bool NextDirection(const LineSearchMinimizer::State& previous, + const LineSearchMinimizer::State& current, + Vector* search_direction) { + low_rank_inverse_hessian_.Update( + previous.search_direction * previous.step_size, + current.gradient - previous.gradient); + search_direction->setZero(); + low_rank_inverse_hessian_.RightMultiply(current.gradient.data(), + search_direction->data()); + *search_direction *= -1.0; + return true; + } + + private: + LowRankInverseHessian low_rank_inverse_hessian_; +}; + +LineSearchDirection* +LineSearchDirection::Create(LineSearchDirection::Options& options) { + if (options.type == STEEPEST_DESCENT) { + return new SteepestDescent; + } + + if (options.type == NONLINEAR_CONJUGATE_GRADIENT) { + return new NonlinearConjugateGradient( + options.nonlinear_conjugate_gradient_type, + options.function_tolerance); + } + + if (options.type == ceres::LBFGS) { + return new ceres::internal::LBFGS(options.num_parameters, + options.max_lbfgs_rank); + } + + LOG(ERROR) << "Unknown line search direction type: " << options.type; + return NULL; +} + +} // namespace internal +} // namespace ceres
diff --git a/internal/ceres/line_search_direction.h b/internal/ceres/line_search_direction.h new file mode 100644 index 0000000..d92c362 --- /dev/null +++ b/internal/ceres/line_search_direction.h
@@ -0,0 +1,71 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2012 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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: sameeragarwal@google.com (Sameer Agarwal) + +#ifndef CERES_INTERNAL_LINE_SEARCH_DIRECTION_H_ +#define CERES_INTERNAL_LINE_SEARCH_DIRECTION_H_ + +#include "ceres/internal/eigen.h" +#include "ceres/line_search_minimizer.h" +#include "ceres/types.h" + +namespace ceres { +namespace internal { + +class LineSearchDirection { + public: + struct Options { + Options() + : num_parameters(0), + type(LBFGS), + nonlinear_conjugate_gradient_type(FLETCHER_REEVES), + function_tolerance(1e-12), + max_lbfgs_rank(20) { + } + + int num_parameters; + LineSearchDirectionType type; + NonlinearConjugateGradientType nonlinear_conjugate_gradient_type; + double function_tolerance; + int max_lbfgs_rank; + }; + + static LineSearchDirection* Create(Options& options); + + virtual ~LineSearchDirection() {} + virtual bool NextDirection(const LineSearchMinimizer::State& previous, + const LineSearchMinimizer::State& current, + Vector* search_direction) = 0; + +}; + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_LINE_SEARCH_DIRECTION_H_
diff --git a/internal/ceres/line_search_minimizer.cc b/internal/ceres/line_search_minimizer.cc index a52d4de..3392194 100644 --- a/internal/ceres/line_search_minimizer.cc +++ b/internal/ceres/line_search_minimizer.cc
@@ -43,18 +43,17 @@ #include <algorithm> #include <cstdlib> #include <cmath> -#include <cstring> -#include <limits> #include <string> #include <vector> #include "Eigen/Dense" #include "ceres/array_utils.h" -#include "ceres/lbfgs.h" #include "ceres/evaluator.h" #include "ceres/internal/eigen.h" +#include "ceres/internal/port.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/line_search.h" +#include "ceres/line_search_direction.h" #include "ceres/stringprintf.h" #include "ceres/types.h" #include "ceres/wall_time.h" @@ -65,35 +64,31 @@ namespace { // Small constant for various floating point issues. const double kEpsilon = 1e-12; -} // namespace - -// Execute the list of IterationCallbacks sequentially. If any one of -// the callbacks does not return SOLVER_CONTINUE, then stop and return -// its status. -CallbackReturnType LineSearchMinimizer::RunCallbacks( - const IterationSummary& iteration_summary) { - for (int i = 0; i < options_.callbacks.size(); ++i) { - const CallbackReturnType status = - (*options_.callbacks[i])(iteration_summary); - if (status != SOLVER_CONTINUE) { - return status; - } +bool Evaluate(Evaluator* evaluator, + const Vector& x, + LineSearchMinimizer::State* state) { + const bool status = evaluator->Evaluate(x.data(), + &(state->cost), + NULL, + state->gradient.data(), + NULL); + if (status) { + state->gradient_squared_norm = state->gradient.squaredNorm(); + state->gradient_max_norm = state->gradient.lpNorm<Eigen::Infinity>(); } - return SOLVER_CONTINUE; + + return status; } -void LineSearchMinimizer::Init(const Minimizer::Options& options) { - options_ = options; -} +} // namespace void LineSearchMinimizer::Minimize(const Minimizer::Options& options, double* parameters, Solver::Summary* summary) { double start_time = WallTimeInSeconds(); double iteration_start_time = start_time; - Init(options); - Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator); + Evaluator* evaluator = CHECK_NOTNULL(options.evaluator); const int num_parameters = evaluator->NumParameters(); const int num_effective_parameters = evaluator->NumEffectiveParameters(); @@ -103,21 +98,12 @@ VectorRef x(parameters, num_parameters); - Vector gradient(num_effective_parameters); - double gradient_squared_norm; - Vector previous_gradient(num_effective_parameters); - Vector gradient_change(num_effective_parameters); - double previous_gradient_squared_norm = 0.0; - - Vector search_direction(num_effective_parameters); - Vector previous_search_direction(num_effective_parameters); + State current_state(num_parameters, num_effective_parameters); + State previous_state(num_parameters, num_effective_parameters); Vector delta(num_effective_parameters); Vector x_plus_delta(num_parameters); - double directional_derivative = 0.0; - double previous_directional_derivative = 0.0; - IterationSummary iteration_summary; iteration_summary.iteration = 0; iteration_summary.step_is_valid = false; @@ -129,31 +115,28 @@ iteration_summary.step_solver_time_in_seconds = 0; // Do initial cost and Jacobian evaluation. - double cost = 0.0; - double previous_cost = 0.0; - if (!evaluator->Evaluate(x.data(), &cost, NULL, gradient.data(), NULL)) { + if (!Evaluate(evaluator, x, ¤t_state)) { LOG(WARNING) << "Terminating: Cost and gradient evaluation failed."; summary->termination_type = NUMERICAL_FAILURE; return; } - gradient_squared_norm = gradient.squaredNorm(); - iteration_summary.cost = cost + summary->fixed_cost; - iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>(); + iteration_summary.cost = current_state.cost + summary->fixed_cost; + iteration_summary.gradient_max_norm = current_state.gradient_max_norm; // The initial gradient max_norm is bounded from below so that we do // not divide by zero. - const double gradient_max_norm_0 = + const double initial_gradient_max_norm = max(iteration_summary.gradient_max_norm, kEpsilon); const double absolute_gradient_tolerance = - options_.gradient_tolerance * gradient_max_norm_0; + options.gradient_tolerance * initial_gradient_max_norm; if (iteration_summary.gradient_max_norm <= absolute_gradient_tolerance) { summary->termination_type = GRADIENT_TOLERANCE; VLOG(1) << "Terminating: Gradient tolerance reached." << "Relative gradient max norm: " - << iteration_summary.gradient_max_norm / gradient_max_norm_0 - << " <= " << options_.gradient_tolerance; + << iteration_summary.gradient_max_norm / initial_gradient_max_norm + << " <= " << options.gradient_tolerance; return; } @@ -164,23 +147,14 @@ + summary->preprocessor_time_in_seconds; summary->iterations.push_back(iteration_summary); - // Call the various callbacks. TODO(sameeragarwal): Here and in - // trust_region_minimizer make this into a function that can be - // shared. - switch (RunCallbacks(iteration_summary)) { - case SOLVER_TERMINATE_SUCCESSFULLY: - summary->termination_type = USER_SUCCESS; - VLOG(1) << "Terminating: User callback returned USER_SUCCESS."; - return; - case SOLVER_ABORT: - summary->termination_type = USER_ABORT; - VLOG(1) << "Terminating: User callback returned USER_ABORT."; - return; - case SOLVER_CONTINUE: - break; - default: - LOG(FATAL) << "Unknown type of user callback status"; - } + LineSearchDirection::Options line_search_direction_options; + line_search_direction_options.num_parameters = num_effective_parameters; + line_search_direction_options.type = options.line_search_direction_type; + line_search_direction_options.nonlinear_conjugate_gradient_type = + options.nonlinear_conjugate_gradient_type; + line_search_direction_options.max_lbfgs_rank = options.max_lbfgs_rank; + scoped_ptr<LineSearchDirection> line_search_direction( + LineSearchDirection::Create(line_search_direction_options)); LineSearchFunction line_search_function(evaluator); LineSearch::Options line_search_options; @@ -191,14 +165,13 @@ ArmijoLineSearch line_search; LineSearch::Summary line_search_summary; - scoped_ptr<LBFGS> lbfgs; - if (options_.line_search_direction_type == ceres::LBFGS) { - lbfgs.reset(new LBFGS(num_effective_parameters, options_.max_lbfgs_rank)); - } - while (true) { + if (!RunCallbacks(options.callbacks, iteration_summary, summary)) { + return; + } + iteration_start_time = WallTimeInSeconds(); - if (iteration_summary.iteration >= options_.max_num_iterations) { + if (iteration_summary.iteration >= options.max_num_iterations) { summary->termination_type = NO_CONVERGENCE; VLOG(1) << "Terminating: Maximum number of iterations reached."; break; @@ -206,175 +179,98 @@ const double total_solver_time = iteration_start_time - start_time + summary->preprocessor_time_in_seconds; - if (total_solver_time >= options_.max_solver_time_in_seconds) { + if (total_solver_time >= options.max_solver_time_in_seconds) { summary->termination_type = NO_CONVERGENCE; VLOG(1) << "Terminating: Maximum solver time reached."; break; } - previous_search_direction = search_direction; - iteration_summary = IterationSummary(); iteration_summary.iteration = summary->iterations.back().iteration + 1; - iteration_summary.step_is_valid = false; - iteration_summary.step_is_successful = false; + bool line_search_status = true; if (iteration_summary.iteration == 1) { - search_direction = -gradient; - directional_derivative = -gradient_squared_norm; + current_state.search_direction = -current_state.gradient; } else { - if (lbfgs.get() != NULL) { - lbfgs->Update(delta, gradient_change); - } - - // TODO(sameeragarwal): This should probably be refactored into - // a set of functions. But we will do that once things settle - // down in this solver. - switch (options_.line_search_direction_type) { - case STEEPEST_DESCENT: - search_direction = -gradient; - directional_derivative = -gradient_squared_norm; - break; - - case NONLINEAR_CONJUGATE_GRADIENT: - { - double beta = 0.0; - - switch (options_.nonlinear_conjugate_gradient_type) { - case FLETCHER_REEVES: - beta = gradient.squaredNorm() / - previous_gradient_squared_norm; - break; - - case POLAK_RIBIRERE: - gradient_change = gradient - previous_gradient; - beta = gradient.dot(gradient_change) / - previous_gradient_squared_norm; - break; - - case HESTENES_STIEFEL: - gradient_change = gradient - previous_gradient; - beta = gradient.dot(gradient_change) / - previous_search_direction.dot(gradient_change); - break; - - default: - LOG(FATAL) << "Unknown nonlinear conjugate gradient type: " - << options_.nonlinear_conjugate_gradient_type; - } - - search_direction = -gradient + beta * previous_search_direction; - } - - directional_derivative = gradient.dot(search_direction); - if (directional_derivative > -options.function_tolerance) { - LOG(WARNING) << "Restarting non-linear conjugate gradients: " - << directional_derivative; - search_direction = -gradient; - directional_derivative = -gradient_squared_norm; - } - break; - - case ceres::LBFGS: - search_direction.setZero(); - lbfgs->RightMultiply(gradient.data(), search_direction.data()); - search_direction *= -1.0; - directional_derivative = gradient.dot(search_direction); - break; - - default: - LOG(FATAL) << "Unknown line search direction type: " - << options_.line_search_direction_type; - } + line_search_status = line_search_direction->NextDirection( + previous_state, + current_state, + ¤t_state.search_direction); } + if (!line_search_status) { + LOG(WARNING) << "Line search direction computation failed. " + "Resorting to steepest descent."; + current_state.search_direction = -current_state.gradient; + } + + line_search_function.Init(x, current_state.search_direction); + current_state.directional_derivative = + current_state.gradient.dot(current_state.search_direction); + // TODO(sameeragarwal): Refactor this into its own object and add // explanations for the various choices. const double initial_step_size = (iteration_summary.iteration == 1) - ? min(1.0, 1.0 / gradient.lpNorm<Eigen::Infinity>()) - : min(1.0, 2.0 * (cost - previous_cost) / directional_derivative); + ? min(1.0, 1.0 / current_state.gradient_max_norm) + : min(1.0, 2.0 * (current_state.cost - previous_state.cost) / + current_state.directional_derivative); - previous_cost = cost; - previous_gradient = gradient; - previous_gradient_squared_norm = gradient_squared_norm; - previous_directional_derivative = directional_derivative; - - line_search_function.Init(x, search_direction); line_search.Search(line_search_options, initial_step_size, - cost, - directional_derivative, + current_state.cost, + current_state.directional_derivative, &line_search_summary); - delta = line_search_summary.optimal_step_size * search_direction; + current_state.step_size = line_search_summary.optimal_step_size; + delta = current_state.step_size * current_state.search_direction; + + previous_state = current_state; + // TODO(sameeragarwal): Collect stats. if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data()) || - !evaluator->Evaluate(x_plus_delta.data(), - &cost, - NULL, - gradient.data(), - NULL)) { + !Evaluate(evaluator, x_plus_delta, ¤t_state)) { LOG(WARNING) << "Evaluation failed."; - cost = previous_cost; - gradient = previous_gradient; } else { x = x_plus_delta; - gradient_squared_norm = gradient.squaredNorm(); } - - iteration_summary.cost = cost + summary->fixed_cost; - iteration_summary.cost_change = previous_cost - cost; - iteration_summary.step_norm = delta.norm(); - iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>(); - iteration_summary.step_is_valid = true; - iteration_summary.step_is_successful = true; - iteration_summary.step_norm = delta.norm(); - iteration_summary.step_size = line_search_summary.optimal_step_size; - iteration_summary.line_search_function_evaluations = - line_search_summary.num_evaluations; - + iteration_summary.gradient_max_norm = current_state.gradient_max_norm; if (iteration_summary.gradient_max_norm <= absolute_gradient_tolerance) { summary->termination_type = GRADIENT_TOLERANCE; VLOG(1) << "Terminating: Gradient tolerance reached." << "Relative gradient max norm: " - << iteration_summary.gradient_max_norm / gradient_max_norm_0 - << " <= " << options_.gradient_tolerance; + << iteration_summary.gradient_max_norm / initial_gradient_max_norm + << " <= " << options.gradient_tolerance; break; } + iteration_summary.cost_change = previous_state.cost - current_state.cost; const double absolute_function_tolerance = - options_.function_tolerance * previous_cost; + options.function_tolerance * previous_state.cost; if (fabs(iteration_summary.cost_change) < absolute_function_tolerance) { VLOG(1) << "Terminating. Function tolerance reached. " << "|cost_change|/cost: " - << fabs(iteration_summary.cost_change) / previous_cost - << " <= " << options_.function_tolerance; + << fabs(iteration_summary.cost_change) / previous_state.cost + << " <= " << options.function_tolerance; summary->termination_type = FUNCTION_TOLERANCE; return; } + iteration_summary.cost = current_state.cost + summary->fixed_cost; + iteration_summary.step_norm = delta.norm(); + iteration_summary.step_is_valid = true; + iteration_summary.step_is_successful = true; + iteration_summary.step_norm = delta.norm(); + iteration_summary.step_size = current_state.step_size; + iteration_summary.line_search_function_evaluations = + line_search_summary.num_evaluations; iteration_summary.iteration_time_in_seconds = WallTimeInSeconds() - iteration_start_time; iteration_summary.cumulative_time_in_seconds = WallTimeInSeconds() - start_time + summary->preprocessor_time_in_seconds; - summary->iterations.push_back(iteration_summary); - switch (RunCallbacks(iteration_summary)) { - case SOLVER_TERMINATE_SUCCESSFULLY: - summary->termination_type = USER_SUCCESS; - VLOG(1) << "Terminating: User callback returned USER_SUCCESS."; - return; - case SOLVER_ABORT: - summary->termination_type = USER_ABORT; - VLOG(1) << "Terminating: User callback returned USER_ABORT."; - return; - case SOLVER_CONTINUE: - break; - default: - LOG(FATAL) << "Unknown type of user callback status"; - } + summary->iterations.push_back(iteration_summary); } }
diff --git a/internal/ceres/line_search_minimizer.h b/internal/ceres/line_search_minimizer.h index a9556fd..f82f139 100644 --- a/internal/ceres/line_search_minimizer.h +++ b/internal/ceres/line_search_minimizer.h
@@ -34,6 +34,8 @@ #include "ceres/minimizer.h" #include "ceres/solver.h" #include "ceres/types.h" +#include "ceres/internal/eigen.h" +#include "glog/logging.h" namespace ceres { namespace internal { @@ -43,17 +45,33 @@ // For example usage, see SolverImpl::Minimize. class LineSearchMinimizer : public Minimizer { public: + struct State { + State(int num_parameters, + int num_effective_parameters) + : cost(0.0), + gradient(num_effective_parameters), + gradient_squared_norm(0.0), + search_direction(num_effective_parameters), + directional_derivative(0.0), + step_size(0.0) { + } + + double cost; + Vector gradient; + double gradient_squared_norm; + double gradient_max_norm; + Vector search_direction; + double directional_derivative; + double step_size; + }; + ~LineSearchMinimizer() {} virtual void Minimize(const Minimizer::Options& options, double* parameters, Solver::Summary* summary); - - private: - void Init(const Minimizer::Options& options); - CallbackReturnType RunCallbacks(const IterationSummary& iteration_summary); - Minimizer::Options options_; }; } // namespace internal } // namespace ceres + #endif // CERES_INTERNAL_LINE_SEARCH_MINIMIZER_H_
diff --git a/internal/ceres/lbfgs.cc b/internal/ceres/low_rank_inverse_hessian.cc similarity index 90% rename from internal/ceres/lbfgs.cc rename to internal/ceres/low_rank_inverse_hessian.cc index 48067c8..3fe113f 100644 --- a/internal/ceres/lbfgs.cc +++ b/internal/ceres/low_rank_inverse_hessian.cc
@@ -28,14 +28,15 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) -#include "glog/logging.h" -#include "ceres/lbfgs.h" #include "ceres/internal/eigen.h" +#include "ceres/low_rank_inverse_hessian.h" +#include "glog/logging.h" namespace ceres { namespace internal { -LBFGS::LBFGS(int num_parameters, int max_num_corrections) +LowRankInverseHessian::LowRankInverseHessian(int num_parameters, + int max_num_corrections) : num_parameters_(num_parameters), max_num_corrections_(max_num_corrections), num_corrections_(0), @@ -45,7 +46,8 @@ delta_x_dot_delta_gradient_(max_num_corrections) { } -bool LBFGS::Update(const Vector& delta_x, const Vector& delta_gradient) { +bool LowRankInverseHessian::Update(const Vector& delta_x, + const Vector& delta_gradient) { const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient); if (delta_x_dot_delta_gradient <= 1e-10) { VLOG(2) << "Skipping LBFGS Update. " << delta_x_dot_delta_gradient; @@ -79,7 +81,8 @@ return true; } -void LBFGS::RightMultiply(const double* x_ptr, double* y_ptr) const { +void LowRankInverseHessian::RightMultiply(const double* x_ptr, + double* y_ptr) const { ConstVectorRef gradient(x_ptr, num_parameters_); VectorRef search_direction(y_ptr, num_parameters_);
diff --git a/internal/ceres/lbfgs.h b/internal/ceres/low_rank_inverse_hessian.h similarity index 71% rename from internal/ceres/lbfgs.h rename to internal/ceres/low_rank_inverse_hessian.h index f6281d8..5d59c14 100644 --- a/internal/ceres/lbfgs.h +++ b/internal/ceres/low_rank_inverse_hessian.h
@@ -29,14 +29,10 @@ // Author: sameeragarwal@google.com (Sameer Agarwal) // // Limited memory positive definite approximation to the inverse -// Hessian, using the work of -// -// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with Limited -// Storage". Mathematics of Computation 35 (151): 773–782. -// -// Byrd, R. H.; Nocedal, J.; Schnabel, R. B. (1994). -// "Representations of Quasi-Newton Matrices and their use in -// Limited Memory Methods". Mathematical Programming 63 (4): +// Hessian, using the LBFGS algorithm + +#ifndef CERES_INTERNAL_LOW_RANK_INVERSE_HESSIAN_H_ +#define CERES_INTERNAL_LOW_RANK_INVERSE_HESSIAN_H_ #include "ceres/internal/eigen.h" #include "ceres/linear_operator.h" @@ -44,22 +40,39 @@ namespace ceres { namespace internal { -// Right multiplying with an LBFGS operator is equivalent to -// multiplying by the inverse Hessian. -class LBFGS : public LinearOperator { +// LowRankInverseHessian is a positive definite approximation to the +// Hessian using the limited memory variant of the +// Broyden-Fletcher-Goldfarb-Shanno (BFGS)secant formula for +// approximating the Hessian. +// +// Other update rules like the Davidon-Fletcher-Powell (DFP) are +// possible, but the BFGS rule is considered the best performing one. +// +// The limited memory variant was developed by Nocedal and further +// enhanced with scaling rule by Byrd, Nocedal and Schanbel. +// +// Nocedal, J. (1980). "Updating Quasi-Newton Matrices with Limited +// Storage". Mathematics of Computation 35 (151): 773–782. +// +// Byrd, R. H.; Nocedal, J.; Schnabel, R. B. (1994). +// "Representations of Quasi-Newton Matrices and their use in +// Limited Memory Methods". Mathematical Programming 63 (4): +class LowRankInverseHessian : public LinearOperator { public: // num_parameters is the row/column size of the Hessian. // max_num_corrections is the rank of the Hessian approximation. // The approximation uses: // 2 * max_num_corrections * num_parameters + max_num_corrections // doubles. - LBFGS(int num_parameters, int max_num_corrections); - virtual ~LBFGS() {} + LowRankInverseHessian(int num_parameters, int max_num_corrections); + virtual ~LowRankInverseHessian() {} - // Update the low rank approximation, i.e. store delta_x and - // delta_gradient, and get rid of the oldest delta_x and - // delta_gradient vectors if the number of corrections is already - // equal to max_num_corrections. + // Update the low rank approximation. delta_x is the change in the + // domain of Hessian, and delta_gradient is the change in the + // gradient. The update copies the delta_x and delta_gradient + // vectors, and gets rid of the oldest delta_x and delta_gradient + // vectors if the number of corrections is already equal to + // max_num_corrections. bool Update(const Vector& delta_x, const Vector& delta_gradient); // LinearOperator interface @@ -82,3 +95,5 @@ } // namespace internal } // namespace ceres + +#endif // CERES_INTERNAL_LOW_RANK_INVERSE_HESSIAN_H_
diff --git a/internal/ceres/minimizer.cc b/internal/ceres/minimizer.cc new file mode 100644 index 0000000..912ac99 --- /dev/null +++ b/internal/ceres/minimizer.cc
@@ -0,0 +1,67 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2010, 2011, 2012 Google Inc. All rights reserved. +// http://code.google.com/p/ceres-solver/ +// +// 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: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/minimizer.h" +#include "ceres/types.h" +#include "glog/logging.h" + +namespace ceres { +namespace internal { + +Minimizer::~Minimizer() {}; + +bool Minimizer::RunCallbacks(const vector<IterationCallback*> callbacks, + const IterationSummary& iteration_summary, + Solver::Summary* summary) { + CallbackReturnType status = SOLVER_CONTINUE; + int i = 0; + while (status == SOLVER_CONTINUE && i < callbacks.size()) { + status = (*callbacks[i])(iteration_summary); + ++i; + } + switch (status) { + case SOLVER_CONTINUE: + return true; + case SOLVER_TERMINATE_SUCCESSFULLY: + summary->termination_type = USER_SUCCESS; + VLOG(1) << "Terminating: User callback returned USER_SUCCESS."; + return false; + case SOLVER_ABORT: + summary->termination_type = USER_ABORT; + VLOG(1) << "Terminating: User callback returned USER_ABORT."; + return false; + default: + LOG(FATAL) << "Unknown type of user callback status"; + } + return false; +} + +} // namespace internal +} // namespace ceres
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h index 2be3a3b..708974d 100644 --- a/internal/ceres/minimizer.h +++ b/internal/ceres/minimizer.h
@@ -143,8 +143,11 @@ Minimizer* inner_iteration_minimizer; }; - virtual ~Minimizer() {} + static bool RunCallbacks(const vector<IterationCallback*> callbacks, + const IterationSummary& iteration_summary, + Solver::Summary* summary); + virtual ~Minimizer(); // Note: The minimizer is expected to update the state of the // parameters array every iteration. This is required for the // StateUpdatingCallback to work.
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc index db2641d..5ea9374 100644 --- a/internal/ceres/trust_region_minimizer.cc +++ b/internal/ceres/trust_region_minimizer.cc
@@ -58,21 +58,6 @@ const double kEpsilon = 1e-12; } // namespace -// Execute the list of IterationCallbacks sequentially. If any one of -// the callbacks does not return SOLVER_CONTINUE, then stop and return -// its status. -CallbackReturnType TrustRegionMinimizer::RunCallbacks( - const IterationSummary& iteration_summary) { - for (int i = 0; i < options_.callbacks.size(); ++i) { - const CallbackReturnType status = - (*options_.callbacks[i])(iteration_summary); - if (status != SOLVER_CONTINUE) { - return status; - } - } - return SOLVER_CONTINUE; -} - // Compute a scaling vector that is used to improve the conditioning // of the Jacobian. void TrustRegionMinimizer::EstimateScale(const SparseMatrix& jacobian, @@ -182,16 +167,16 @@ // The initial gradient max_norm is bounded from below so that we do // not divide by zero. - const double gradient_max_norm_0 = + const double initial_gradient_max_norm = max(iteration_summary.gradient_max_norm, kEpsilon); const double absolute_gradient_tolerance = - options_.gradient_tolerance * gradient_max_norm_0; + options_.gradient_tolerance * initial_gradient_max_norm; if (iteration_summary.gradient_max_norm <= absolute_gradient_tolerance) { summary->termination_type = GRADIENT_TOLERANCE; VLOG(1) << "Terminating: Gradient tolerance reached." << "Relative gradient max norm: " - << iteration_summary.gradient_max_norm / gradient_max_norm_0 + << iteration_summary.gradient_max_norm / initial_gradient_max_norm << " <= " << options_.gradient_tolerance; return; } @@ -203,24 +188,12 @@ + summary->preprocessor_time_in_seconds; summary->iterations.push_back(iteration_summary); - // Call the various callbacks. - switch (RunCallbacks(iteration_summary)) { - case SOLVER_TERMINATE_SUCCESSFULLY: - summary->termination_type = USER_SUCCESS; - VLOG(1) << "Terminating: User callback returned USER_SUCCESS."; - return; - case SOLVER_ABORT: - summary->termination_type = USER_ABORT; - VLOG(1) << "Terminating: User callback returned USER_ABORT."; - return; - case SOLVER_CONTINUE: - break; - default: - LOG(FATAL) << "Unknown type of user callback status"; - } - - int num_consecutive_invalid_steps = 0; + int num_consecutive_invalid_steps = 0; while (true) { + if (!RunCallbacks(options.callbacks, iteration_summary, summary)) { + return; + } + iteration_start_time = WallTimeInSeconds(); if (iteration_summary.iteration >= options_.max_num_iterations) { summary->termination_type = NO_CONVERGENCE; @@ -356,7 +329,7 @@ new_cost = x_plus_delta_cost; } else { x_plus_delta = inner_iteration_x; - // Bost the model_cost_change, since the inner iteration + // Boost the model_cost_change, since the inner iteration // improvements are not accounted for by the trust region. model_cost_change += x_plus_delta_cost - new_cost; VLOG(2) << "Inner iteration succeeded; current cost: " << cost @@ -461,7 +434,7 @@ summary->termination_type = GRADIENT_TOLERANCE; VLOG(1) << "Terminating: Gradient tolerance reached." << "Relative gradient max norm: " - << iteration_summary.gradient_max_norm / gradient_max_norm_0 + << iteration_summary.gradient_max_norm / initial_gradient_max_norm << " <= " << options_.gradient_tolerance; return; } @@ -534,21 +507,6 @@ WallTimeInSeconds() - start_time + summary->preprocessor_time_in_seconds; summary->iterations.push_back(iteration_summary); - - switch (RunCallbacks(iteration_summary)) { - case SOLVER_TERMINATE_SUCCESSFULLY: - summary->termination_type = USER_SUCCESS; - VLOG(1) << "Terminating: User callback returned USER_SUCCESS."; - return; - case SOLVER_ABORT: - summary->termination_type = USER_ABORT; - VLOG(1) << "Terminating: User callback returned USER_ABORT."; - return; - case SOLVER_CONTINUE: - break; - default: - LOG(FATAL) << "Unknown type of user callback status"; - } } }
diff --git a/internal/ceres/trust_region_minimizer.h b/internal/ceres/trust_region_minimizer.h index 49f5733..173b0f3 100644 --- a/internal/ceres/trust_region_minimizer.h +++ b/internal/ceres/trust_region_minimizer.h
@@ -52,7 +52,6 @@ private: void Init(const Minimizer::Options& options); void EstimateScale(const SparseMatrix& jacobian, double* scale) const; - CallbackReturnType RunCallbacks(const IterationSummary& iteration_summary); bool MaybeDumpLinearLeastSquaresProblem( const int iteration, const SparseMatrix* jacobian, const double* residuals,
diff --git a/jni/Android.mk b/jni/Android.mk index 349a9c1..f705740 100644 --- a/jni/Android.mk +++ b/jni/Android.mk
@@ -123,13 +123,16 @@ $(CERES_SRC_PATH)/iterative_schur_complement_solver.cc \ $(CERES_SRC_PATH)/levenberg_marquardt_strategy.cc \ $(CERES_SRC_PATH)/line_search.cc \ + $(CERES_SRC_PATH)/line_search_direction.cc \ $(CERES_SRC_PATH)/line_search_minimizer.cc \ $(CERES_SRC_PATH)/linear_least_squares_problems.cc \ $(CERES_SRC_PATH)/linear_operator.cc \ $(CERES_SRC_PATH)/linear_solver.cc \ $(CERES_SRC_PATH)/local_parameterization.cc \ $(CERES_SRC_PATH)/loss_function.cc \ + $(CERES_SRC_PATH)/low_rank_inverse_hessian.cc \ $(CERES_SRC_PATH)/miniglog/glog/logging.cc \ + $(CERES_SRC_PATH)/minimizer.cc \ $(CERES_SRC_PATH)/normal_prior.cc \ $(CERES_SRC_PATH)/parameter_block_ordering.cc \ $(CERES_SRC_PATH)/partitioned_matrix_view.cc \