New Trust region loop.
1. New TrustRegionMinimizer and basic tests for it.
2. New TrustRegionStrategy interface.
3. New LevenbergMarquardtStrategy and tests for it.
4. Updates to SolverImpl to reflect this.
5. Changes to Solver::Options and IterationSummary related to this.
6. Deleted levenberg_marquardt.cc/h/_test.cc
Change-Id: I6c1d1a7c774f014856f9f26263a830aa886e1400
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index 0331f32..718fde9 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -64,15 +64,12 @@
#include "ceres/ceres.h"
DEFINE_string(input, "", "Input File name");
-
DEFINE_string(solver_type, "sparse_schur", "Options are: "
"sparse_schur, dense_schur, iterative_schur, cholesky, "
"dense_qr, and conjugate_gradients");
-
DEFINE_string(preconditioner_type, "jacobi", "Options are: "
"identity, jacobi, schur_jacobi, cluster_jacobi, "
"cluster_tridiagonal");
-
DEFINE_int32(num_iterations, 5, "Number of iterations");
DEFINE_int32(num_threads, 1, "Number of threads");
DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the "
@@ -201,9 +198,9 @@
void SetSolverOptionsFromFlags(BALProblem* bal_problem,
Solver::Options* options) {
+ SetMinimizerOptions(options);
SetLinearSolver(options);
SetOrdering(bal_problem, options);
- SetMinimizerOptions(options);
}
void BuildProblem(BALProblem* bal_problem, Problem* problem) {
diff --git a/include/ceres/iteration_callback.h b/include/ceres/iteration_callback.h
index 88da992..6d7720c 100644
--- a/include/ceres/iteration_callback.h
+++ b/include/ceres/iteration_callback.h
@@ -28,8 +28,9 @@
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
//
-// When an iteration callback is specified, Ceres calls the callback after each
-// optimizer step and pass it an IterationSummary object, defined below.
+// When an iteration callback is specified, Ceres calls the callback
+// after each minimizer step (if the minimizer has not converged) and
+// passes it an IterationSummary object, defined below.
#ifndef CERES_PUBLIC_ITERATION_CALLBACK_H_
#define CERES_PUBLIC_ITERATION_CALLBACK_H_
@@ -44,7 +45,15 @@
// Current iteration number.
int32 iteration;
+ // Step was numerically valid, i.e., all values are finite and the
+ // step reduces the value of the linearized model.
+ //
+ // Note: step_is_valid is false when iteration = 0.
+ bool step_is_valid;
+
// Whether or not the algorithm made progress in this iteration.
+ //
+ // Note: step_is_successful is false when iteration = 0.
bool step_is_successful;
// Value of the objective function.
@@ -66,9 +75,10 @@
// cost and the change in the cost of the linearized approximation.
double relative_decrease;
- // Value of the regularization parameter for Levenberg-Marquardt
- // algorithm at the end of the current iteration.
- double mu;
+ // Size of the trust region at the end of the current iteration. For
+ // the Levenberg-Marquardt algorithm, the regularization parameter
+ // mu = 1.0 / trust_region_radius.
+ double trust_region_radius;
// For the inexact step Levenberg-Marquardt algorithm, this is the
// relative accuracy with which the Newton(LM) step is solved. This
@@ -81,13 +91,15 @@
// Newton step.
int linear_solver_iterations;
- // TODO(sameeragarwal): Change to use a higher precision timer using
- // clock_gettime.
- // Time (in seconds) spent inside the linear least squares solver.
+ // TODO(sameeragarwal): Change the following two to use a higher
+ // precision timer using clock_gettime.
+ //
+ // Time (in seconds) spent inside the minimizer loop in the current
+ // iteration.
int iteration_time_sec;
- // Time (in seconds) spent inside the linear least squares solver.
- int linear_solver_time_sec;
+ // Time (in seconds) spent inside the trust region step solver.
+ int step_solver_time_sec;
};
// Interface for specifying callbacks that are executed at the end of
@@ -133,7 +145,7 @@
// summary.gradient_max_norm,
// summary.step_norm,
// summary.relative_decrease,
-// summary.mu,
+// summary.trust_region_radius,
// summary.eta,
// summary.linear_solver_iterations);
// if (log_to_stdout_) {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 701f0b2..ed4c9b8 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -57,13 +57,17 @@
struct Options {
// Default constructor that sets up a generic sparse problem.
Options() {
- minimizer_type = LEVENBERG_MARQUARDT;
+ trust_region_strategy_type = LEVENBERG_MARQUARDT;
max_num_iterations = 50;
- max_solver_time_sec = 1.0e9;
+ max_solver_time_sec = 1e9;
num_threads = 1;
- tau = 1e-4;
- min_mu = 1e-20;
+ initial_trust_region_radius = 1e4;
+ max_trust_region_radius = 1e16;
+ min_trust_region_radius = 1e-32;
min_relative_decrease = 1e-3;
+ lm_min_diagonal = 1e-6;
+ lm_max_diagonal = 1e32;
+ max_num_consecutive_invalid_steps = 5;
function_tolerance = 1e-6;
gradient_tolerance = 1e-10;
parameter_tolerance = 1e-8;
@@ -93,7 +97,6 @@
return_final_residuals = false;
lsqp_dump_directory = "/tmp";
lsqp_dump_format_type = TEXTFILE;
- crash_and_dump_lsqp_on_failure = false;
check_gradients = false;
gradient_check_relative_precision = 1e-8;
numeric_derivative_relative_step_size = 1e-6;
@@ -102,7 +105,7 @@
// Minimizer options ----------------------------------------
- MinimizerType minimizer_type;
+ TrustRegionStrategyType trust_region_strategy_type;
// Maximum number of iterations for the minimizer to run for.
int max_num_iterations;
@@ -114,27 +117,35 @@
// jacobians.
int num_threads;
- // For Levenberg-Marquardt, the initial value for the
- // regularizer. This is the inversely related to the size of the
- // initial trust region.
- double tau;
+ // Trust region minimizer settings.
+ double initial_trust_region_radius;
+ double max_trust_region_radius;
- // For Levenberg-Marquardt, the minimum value of the
- // regularizer. For well constrained problems there shold be no
- // need to modify the default value, but in some cases, going
- // below a certain minimum reliably triggers rank deficiency in
- // the normal equations. In such cases, the LM solver can
- // oscillate between lowering the value of mu, seeing a numerical
- // failure, and then increasing it making some progress and then
- // reducing it again.
- //
- // In such cases, it is useful to set a higher value for min_mu.
- double min_mu;
+ // Minimizer terminates when the trust region radius becomes
+ // smaller than this value.
+ double min_trust_region_radius;
- // For trust region methods, this is lower threshold for the
- // relative decrease before a step is accepted.
+ // Lower bound for the relative decrease before a step is
+ // accepted.
double min_relative_decrease;
+ // For the Levenberg-Marquadt algorithm, the scaled diagonal of
+ // the normal equations J'J is used to control the size of the
+ // trust region. Extremely small and large values along the
+ // diagonal can make this regularization scheme
+ // fail. lm_max_diagonal and lm_min_diagonal, clamp the values of
+ // diag(J'J) from above and below. In the normal course of
+ // operation, the user should not have to modify these parameters.
+ double lm_min_diagonal;
+ double lm_max_diagonal;
+
+ // Sometimes due to numerical conditioning problems or linear
+ // solver flakiness, the trust region strategy may return a
+ // numerically invalid step that can be fixed by reducing the
+ // trust region size. So the TrustRegionMinimizer allows for a few
+ // successive invalid steps before it declares NUMERICAL_FAILURE.
+ int max_num_consecutive_invalid_steps;
+
// Minimizer terminates when
//
// (new_cost - old_cost) < function_tolerance * old_cost;
@@ -243,15 +254,6 @@
string lsqp_dump_directory;
DumpFormatType lsqp_dump_format_type;
- // Dump the linear least squares problem to disk if the minimizer
- // fails due to NUMERICAL_FAILURE and crash the process. This flag
- // is useful for generating debugging information. The problem is
- // dumped in a file whose name is determined by
- // Solver::Options::lsqp_dump_format.
- //
- // Note: This requires a version of Ceres built with protocol buffers.
- bool crash_and_dump_lsqp_on_failure;
-
// Finite differences options ----------------------------------------------
// Check all jacobians computed by each residual block with finite
@@ -299,10 +301,15 @@
bool update_state_every_iteration;
// Callbacks that are executed at the end of each iteration of the
- // Minimizer. They are executed in the order that they are
- // specified in this vector. By default, parameter blocks are
- // updated only at the end of the optimization, i.e when the
- // Minimizer terminates. This behaviour is controlled by
+ // Minimizer. An iteration may terminate midway, either due to
+ // numerical failures or because one of the convergence tests has
+ // been satisfied. In this case none of the callbacks are
+ // executed.
+
+ // Callbacks are executed in the order that they are specified in
+ // this vector. By default, parameter blocks are updated only at
+ // the end of the optimization, i.e when the Minimizer
+ // terminates. This behaviour is controlled by
// update_state_every_variable. If the user wishes to have access
// to the update parameter blocks when his/her callbacks are
// executed, then set update_state_every_iteration to true.
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 61a9a0d..705097f 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -158,8 +158,8 @@
PER_MINIMIZER_ITERATION
};
-enum MinimizerType {
- LEVENBERG_MARQUARDT
+enum TrustRegionStrategyType {
+ LEVENBERG_MARQUARDT,
};
enum SolverTerminationType {
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 4cdf055..2a0a6af 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -53,7 +53,7 @@
gradient_checking_cost_function.cc
implicit_schur_complement.cc
iterative_schur_complement_solver.cc
- levenberg_marquardt.cc
+ levenberg_marquardt_strategy.cc
linear_least_squares_problems.cc
linear_operator.cc
linear_solver.cc
@@ -78,6 +78,8 @@
split.cc
stringprintf.cc
suitesparse.cc
+ trust_region_minimizer.cc
+ trust_region_strategy.cc
triplet_sparse_matrix.cc
types.cc
visibility_based_preconditioner.cc
@@ -176,7 +178,7 @@
CERES_TEST(implicit_schur_complement)
CERES_TEST(iterative_schur_complement_solver)
CERES_TEST(jet)
- CERES_TEST(levenberg_marquardt)
+ CERES_TEST(levenberg_marquardt_strategy)
CERES_TEST(local_parameterization)
CERES_TEST(loss_function)
CERES_TEST(normal_prior)
@@ -196,6 +198,7 @@
CERES_TEST(solver_impl)
CERES_TEST(symmetric_linear_solver)
CERES_TEST(triplet_sparse_matrix)
+ CERES_TEST(trust_region_minimizer)
CERES_TEST(unsymmetric_linear_solver)
CERES_TEST(visibility)
IF (GFLAGS)
diff --git a/internal/ceres/levenberg_marquardt.cc b/internal/ceres/levenberg_marquardt.cc
deleted file mode 100644
index 1eae624..0000000
--- a/internal/ceres/levenberg_marquardt.cc
+++ /dev/null
@@ -1,588 +0,0 @@
-// 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)
-//
-// Implementation of a simple LM algorithm which uses the step sizing
-// rule of "Methods for Nonlinear Least Squares" by K. Madsen,
-// H.B. Nielsen and O. Tingleff. Available to download from
-//
-// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
-//
-// The basic algorithm described in this note is an exact step
-// algorithm that depends on the Newton(LM) step being solved exactly
-// in each iteration. When a suitable iterative solver is available to
-// solve the Newton(LM) step, the algorithm will automatically switch
-// to an inexact step solution method. This trades some slowdown in
-// convergence for significant savings in solve time and memory
-// usage. Our implementation of the Truncated Newton algorithm follows
-// the discussion and recommendataions in "Stephen G. Nash, A Survey
-// of Truncated Newton Methods, Journal of Computational and Applied
-// Mathematics, 124(1-2), 45-59, 2000.
-
-#include "ceres/levenberg_marquardt.h"
-
-#include <algorithm>
-#include <cstdlib>
-#include <cmath>
-#include <cstring>
-#include <string>
-#include <vector>
-
-#include <glog/logging.h>
-#include "Eigen/Core"
-#include "ceres/array_utils.h"
-#include "ceres/evaluator.h"
-#include "ceres/file.h"
-#include "ceres/linear_least_squares_problems.h"
-#include "ceres/linear_solver.h"
-#include "ceres/matrix_proto.h"
-#include "ceres/sparse_matrix.h"
-#include "ceres/stringprintf.h"
-#include "ceres/internal/eigen.h"
-#include "ceres/internal/scoped_ptr.h"
-#include "ceres/types.h"
-
-namespace ceres {
-namespace internal {
-namespace {
-
-// Numbers for clamping the size of the LM diagonal. The size of these
-// numbers is heuristic. We will probably be adjusting them in the
-// future based on more numerical experience. With jacobi scaling
-// enabled, these numbers should be all but redundant.
-const double kMinLevenbergMarquardtDiagonal = 1e-6;
-const double kMaxLevenbergMarquardtDiagonal = 1e32;
-
-// Small constant for various floating point issues.
-const double kEpsilon = 1e-12;
-
-// Number of times the linear solver should be retried in case of
-// numerical failure. The retries are done by exponentially scaling up
-// mu at each retry. This leads to stronger and stronger
-// regularization making the linear least squares problem better
-// conditioned at each retry.
-const int kMaxLinearSolverRetries = 5;
-
-// D = 1/sqrt(diag(J^T * J))
-void EstimateScale(const SparseMatrix& jacobian, double* D) {
- CHECK_NOTNULL(D);
- jacobian.SquaredColumnNorm(D);
- for (int i = 0; i < jacobian.num_cols(); ++i) {
- D[i] = 1.0 / (kEpsilon + sqrt(D[i]));
- }
-}
-
-// D = diag(J^T * J)
-void LevenbergMarquardtDiagonal(const SparseMatrix& jacobian,
- double* D) {
- CHECK_NOTNULL(D);
- jacobian.SquaredColumnNorm(D);
- for (int i = 0; i < jacobian.num_cols(); ++i) {
- D[i] = min(max(D[i], kMinLevenbergMarquardtDiagonal),
- kMaxLevenbergMarquardtDiagonal);
- }
-}
-
-bool RunCallback(IterationCallback* callback,
- const IterationSummary& iteration_summary,
- Solver::Summary* summary) {
- const CallbackReturnType status = (*callback)(iteration_summary);
- switch (status) {
- case SOLVER_TERMINATE_SUCCESSFULLY:
- summary->termination_type = USER_SUCCESS;
- VLOG(1) << "Terminating on USER_SUCCESS.";
- return false;
- case SOLVER_ABORT:
- summary->termination_type = USER_ABORT;
- VLOG(1) << "Terminating on USER_ABORT.";
- return false;
- case SOLVER_CONTINUE:
- return true;
- default:
- LOG(FATAL) << "Unknown status returned by callback: "
- << status;
- }
-}
-
-} // namespace
-
-LevenbergMarquardt::~LevenbergMarquardt() {}
-
-void LevenbergMarquardt::Minimize(const Minimizer::Options& options,
- Evaluator* evaluator,
- LinearSolver* linear_solver,
- const double* initial_parameters,
- double* final_parameters,
- Solver::Summary* summary) {
- time_t start_time = time(NULL);
- const int num_parameters = evaluator->NumParameters();
- const int num_effective_parameters = evaluator->NumEffectiveParameters();
- const int num_residuals = evaluator->NumResiduals();
-
- summary->termination_type = NO_CONVERGENCE;
- summary->num_successful_steps = 0;
- summary->num_unsuccessful_steps = 0;
-
- // Allocate the various vectors needed by the algorithm.
- memcpy(final_parameters, initial_parameters,
- num_parameters * sizeof(*initial_parameters));
-
- VectorRef x(final_parameters, num_parameters);
- Vector x_new(num_parameters);
-
- Vector lm_step(num_effective_parameters);
- Vector gradient(num_effective_parameters);
- Vector scaled_gradient(num_effective_parameters);
- // Jacobi scaling vector
- Vector scale(num_effective_parameters);
-
- Vector f_model(num_residuals);
- Vector f(num_residuals);
- Vector f_new(num_residuals);
- Vector D(num_parameters);
- Vector muD(num_parameters);
-
- // Ask the Evaluator to create the jacobian matrix. The sparsity
- // pattern of this matrix is going to remain constant, so we only do
- // this once and then re-use this matrix for all subsequent Jacobian
- // computations.
- scoped_ptr<SparseMatrix> jacobian(evaluator->CreateJacobian());
-
- double x_norm = x.norm();
-
- double cost = 0.0;
- D.setOnes();
- f.setZero();
-
- // Do initial cost and Jacobian evaluation.
- if (!evaluator->Evaluate(x.data(), &cost, f.data(), jacobian.get())) {
- LOG(WARNING) << "Failed to compute residuals and Jacobian. "
- << "Terminating.";
- summary->termination_type = NUMERICAL_FAILURE;
- return;
- }
-
- if (options.jacobi_scaling) {
- EstimateScale(*jacobian, scale.data());
- jacobian->ScaleColumns(scale.data());
- } else {
- scale.setOnes();
- }
-
- // This is a poor way to do this computation. Even if fixed_cost is
- // zero, because we are subtracting two possibly large numbers, we
- // are depending on exact cancellation to give us a zero here. But
- // initial_cost and cost have been computed by two different
- // evaluators. One which runs on the whole problem (in
- // solver_impl.cc) in single threaded mode and another which runs
- // here on the reduced problem, so fixed_cost can (and does) contain
- // some numerical garbage with a relative magnitude of 1e-14.
- //
- // The right way to do this, would be to compute the fixed cost on
- // just the set of residual blocks which are held constant and were
- // removed from the original problem when the reduced problem was
- // constructed.
- summary->fixed_cost = summary->initial_cost - cost;
-
- double model_cost = f.squaredNorm() / 2.0;
- double total_cost = summary->fixed_cost + cost;
-
- scaled_gradient.setZero();
- jacobian->LeftMultiply(f.data(), scaled_gradient.data());
- gradient = scaled_gradient.array() / scale.array();
-
- double gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
- // We need the max here to guard againt the gradient being zero.
- const double gradient_max_norm_0 = max(gradient_max_norm, kEpsilon);
- double gradient_tolerance = options.gradient_tolerance * gradient_max_norm_0;
-
- double mu = options.tau;
- double nu = 2.0;
- int iteration = 0;
- double actual_cost_change = 0.0;
- double step_norm = 0.0;
- double relative_decrease = 0.0;
-
- // Insane steps are steps which are not sane, i.e. there is some
- // numerical kookiness going on with them. There are various reasons
- // for this kookiness, some easier to diagnose then others. From the
- // point of view of the non-linear solver, they are steps which
- // cannot be used. We return with NUMERICAL_FAILURE after
- // kMaxLinearSolverRetries consecutive insane steps.
- bool step_is_sane = false;
- int num_consecutive_insane_steps = 0;
-
- // Whether the step resulted in a sufficient decrease in the
- // objective function when compared to the decrease in the value of
- // the lineariztion.
- bool step_is_successful = false;
-
- // Parse the iterations for which to dump the linear problem.
- vector<int> iterations_to_dump = options.lsqp_iterations_to_dump;
- sort(iterations_to_dump.begin(), iterations_to_dump.end());
-
- IterationSummary iteration_summary;
- iteration_summary.iteration = iteration;
- iteration_summary.step_is_successful = false;
- iteration_summary.cost = total_cost;
- iteration_summary.cost_change = actual_cost_change;
- iteration_summary.gradient_max_norm = gradient_max_norm;
- iteration_summary.step_norm = step_norm;
- iteration_summary.relative_decrease = relative_decrease;
- iteration_summary.mu = mu;
- iteration_summary.eta = options.eta;
- iteration_summary.linear_solver_iterations = 0;
- iteration_summary.linear_solver_time_sec = 0.0;
- iteration_summary.iteration_time_sec = (time(NULL) - start_time);
- if (options.logging_type >= PER_MINIMIZER_ITERATION) {
- summary->iterations.push_back(iteration_summary);
- }
-
- // Check if the starting point is an optimum.
- VLOG(2) << "Gradient max norm: " << gradient_max_norm
- << " tolerance: " << gradient_tolerance
- << " ratio: " << gradient_max_norm / gradient_max_norm_0
- << " tolerance: " << options.gradient_tolerance;
- if (gradient_max_norm <= gradient_tolerance) {
- summary->termination_type = GRADIENT_TOLERANCE;
- VLOG(1) << "Terminating on GRADIENT_TOLERANCE. "
- << "Relative gradient max norm: "
- << gradient_max_norm / gradient_max_norm_0
- << " <= " << options.gradient_tolerance;
- return;
- }
-
- // Call the various callbacks.
- for (int i = 0; i < options.callbacks.size(); ++i) {
- if (!RunCallback(options.callbacks[i], iteration_summary, summary)) {
- return;
- }
- }
-
- // We only need the LM diagonal if we are actually going to do at
- // least one iteration of the optimization. So we wait to do it
- // until now.
- LevenbergMarquardtDiagonal(*jacobian, D.data());
-
- while ((iteration < options.max_num_iterations) &&
- (time(NULL) - start_time) <= options.max_solver_time_sec) {
- time_t iteration_start_time = time(NULL);
- step_is_sane = false;
- step_is_successful = false;
-
- IterationSummary iteration_summary;
- // The while loop here is just to provide an easily breakable
- // control structure. We are guaranteed to always exit this loop
- // at the end of one iteration or before.
- while (1) {
- muD = (mu * D).array().sqrt();
- LinearSolver::PerSolveOptions solve_options;
- solve_options.D = muD.data();
- solve_options.q_tolerance = options.eta;
- // Disable r_tolerance checking. Since we only care about
- // termination via the q_tolerance. As Nash and Sofer show,
- // r_tolerance based termination is essentially useless in
- // Truncated Newton methods.
- solve_options.r_tolerance = -1.0;
-
- // Invalidate the output array lm_step, so that we can detect if
- // the linear solver generated numerical garbage. This is known
- // to happen for the DENSE_QR and then DENSE_SCHUR solver when
- // the Jacobin is severly rank deficient and mu is too small.
- InvalidateArray(num_effective_parameters, lm_step.data());
- const time_t linear_solver_start_time = time(NULL);
- LinearSolver::Summary linear_solver_summary =
- linear_solver->Solve(jacobian.get(),
- f.data(),
- solve_options,
- lm_step.data());
- iteration_summary.linear_solver_time_sec =
- (time(NULL) - linear_solver_start_time);
- iteration_summary.linear_solver_iterations =
- linear_solver_summary.num_iterations;
-
- if (binary_search(iterations_to_dump.begin(),
- iterations_to_dump.end(),
- iteration)) {
- CHECK(DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
- iteration,
- options.lsqp_dump_format_type,
- jacobian.get(),
- muD.data(),
- f.data(),
- lm_step.data(),
- options.num_eliminate_blocks))
- << "Tried writing linear least squares problem: "
- << options.lsqp_dump_directory
- << " but failed.";
- }
-
- // We ignore the case where the linear solver did not converge,
- // since the partial solution computed by it still maybe of use,
- // and there is no reason to ignore it, especially since we
- // spent so much time computing it.
- if ((linear_solver_summary.termination_type != TOLERANCE) &&
- (linear_solver_summary.termination_type != MAX_ITERATIONS)) {
- VLOG(1) << "Linear solver failure: retrying with a higher mu";
- break;
- }
-
- if (!IsArrayValid(num_effective_parameters, lm_step.data())) {
- LOG(WARNING) << "Linear solver failure. Failed to compute a finite "
- << "step. Terminating. Please report this to the Ceres "
- << "Solver developers.";
- summary->termination_type = NUMERICAL_FAILURE;
- return;
- }
-
- step_norm = (lm_step.array() * scale.array()).matrix().norm();
-
- // Check step length based convergence. If the step length is
- // too small, then we are done.
- const double step_size_tolerance = options.parameter_tolerance *
- (x_norm + options.parameter_tolerance);
-
- VLOG(2) << "Step size: " << step_norm
- << " tolerance: " << step_size_tolerance
- << " ratio: " << step_norm / step_size_tolerance
- << " tolerance: " << options.parameter_tolerance;
- if (step_norm <= options.parameter_tolerance *
- (x_norm + options.parameter_tolerance)) {
- summary->termination_type = PARAMETER_TOLERANCE;
- VLOG(1) << "Terminating on PARAMETER_TOLERANCE."
- << "Relative step size: " << step_norm / step_size_tolerance
- << " <= " << options.parameter_tolerance;
- return;
- }
-
- Vector delta = -(lm_step.array() * scale.array()).matrix();
- if (!evaluator->Plus(x.data(), delta.data(), x_new.data())) {
- LOG(WARNING) << "Failed to compute Plus(x, delta, x_plus_delta). "
- << "Terminating.";
- summary->termination_type = NUMERICAL_FAILURE;
- return;
- }
-
- double cost_new = 0.0;
- if (!evaluator->Evaluate(x_new.data(), &cost_new, NULL, NULL)) {
- LOG(WARNING) << "Failed to compute the value of the objective "
- << "function. Terminating.";
- summary->termination_type = NUMERICAL_FAILURE;
- return;
- }
-
- f_model.setZero();
- jacobian->RightMultiply(lm_step.data(), f_model.data());
- const double model_cost_new =
- (f.segment(0, num_residuals) - f_model).squaredNorm() / 2;
-
- actual_cost_change = cost - cost_new;
- double model_cost_change = model_cost - model_cost_new;
-
- VLOG(2) << "[Model cost] current: " << model_cost
- << " new : " << model_cost_new
- << " change: " << model_cost_change;
-
- VLOG(2) << "[Nonlinear cost] current: " << cost
- << " new : " << cost_new
- << " change: " << actual_cost_change
- << " relative change: " << fabs(actual_cost_change) / cost
- << " tolerance: " << options.function_tolerance;
-
- // In exact arithmetic model_cost_change should never be
- // negative. But due to numerical precision issues, we may end up
- // with a small negative number. model_cost_change which are
- // negative and large in absolute value are indicative of a
- // numerical failure in the solver.
- if (model_cost_change < -kEpsilon) {
- VLOG(1) << "Model cost change is negative.\n"
- << "Current : " << model_cost
- << " new : " << model_cost_new
- << " change: " << model_cost_change << "\n";
- break;
- }
-
- // If we have reached this far, then we are willing to trust the
- // numerical quality of the step.
- step_is_sane = true;
- num_consecutive_insane_steps = 0;
-
- // Check function value based convergence.
- if (fabs(actual_cost_change) < options.function_tolerance * cost) {
- VLOG(1) << "Termination on FUNCTION_TOLERANCE."
- << " Relative cost change: " << fabs(actual_cost_change) / cost
- << " tolerance: " << options.function_tolerance;
- summary->termination_type = FUNCTION_TOLERANCE;
- return;
- }
-
- // Clamp model_cost_change at kEpsilon from below.
- if (model_cost_change < kEpsilon) {
- VLOG(1) << "Clamping model cost change " << model_cost_change
- << " to " << kEpsilon;
- model_cost_change = kEpsilon;
- }
-
- relative_decrease = actual_cost_change / model_cost_change;
- VLOG(2) << "actual_cost_change / model_cost_change = "
- << relative_decrease;
-
- if (relative_decrease < options.min_relative_decrease) {
- VLOG(2) << "Unsuccessful step.";
- break;
- }
-
- VLOG(2) << "Successful step.";
-
- ++summary->num_successful_steps;
- x = x_new;
- x_norm = x.norm();
-
- if (!evaluator->Evaluate(x.data(), &cost, f.data(), jacobian.get())) {
- LOG(WARNING) << "Failed to compute residuals and jacobian. "
- << "Terminating.";
- summary->termination_type = NUMERICAL_FAILURE;
- return;
- }
-
- if (options.jacobi_scaling) {
- jacobian->ScaleColumns(scale.data());
- }
-
- model_cost = f.squaredNorm() / 2.0;
- LevenbergMarquardtDiagonal(*jacobian, D.data());
- scaled_gradient.setZero();
- jacobian->LeftMultiply(f.data(), scaled_gradient.data());
- gradient = scaled_gradient.array() / scale.array();
- gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
-
- // Check gradient based convergence.
- VLOG(2) << "Gradient max norm: " << gradient_max_norm
- << " tolerance: " << gradient_tolerance
- << " ratio: " << gradient_max_norm / gradient_max_norm_0
- << " tolerance: " << options.gradient_tolerance;
- if (gradient_max_norm <= gradient_tolerance) {
- summary->termination_type = GRADIENT_TOLERANCE;
- VLOG(1) << "Terminating on GRADIENT_TOLERANCE. "
- << "Relative gradient max norm: "
- << gradient_max_norm / gradient_max_norm_0
- << " <= " << options.gradient_tolerance
- << " (tolerance).";
- return;
- }
-
- mu = mu * max(1.0 / 3.0, 1 - pow(2 * relative_decrease - 1, 3));
- mu = std::max(options.min_mu, mu);
- nu = 2.0;
- step_is_successful = true;
- break;
- }
-
- if (!step_is_sane) {
- ++num_consecutive_insane_steps;
- }
-
- if (num_consecutive_insane_steps == kMaxLinearSolverRetries) {
- summary->termination_type = NUMERICAL_FAILURE;
- VLOG(1) << "Too many consecutive retries; ending with numerical fail.";
-
- if (!options.crash_and_dump_lsqp_on_failure) {
- return;
- }
-
- // Dump debugging information to disk.
- CHECK(options.lsqp_dump_format_type == TEXTFILE ||
- options.lsqp_dump_format_type == PROTOBUF)
- << "Dumping the linear least squares problem on crash "
- << "requires Solver::Options::lsqp_dump_format_type to be "
- << "PROTOBUF or TEXTFILE.";
-
- if (DumpLinearLeastSquaresProblem(options.lsqp_dump_directory,
- iteration,
- options.lsqp_dump_format_type,
- jacobian.get(),
- muD.data(),
- f.data(),
- lm_step.data(),
- options.num_eliminate_blocks)) {
- LOG(FATAL) << "Linear least squares problem saved to: "
- << options.lsqp_dump_directory
- << ". Please provide this to the Ceres developers for "
- << " debugging along with the v=2 log.";
- } else {
- LOG(FATAL) << "Tried writing linear least squares problem: "
- << options.lsqp_dump_directory
- << " but failed.";
- }
- }
-
- if (!step_is_successful) {
- // Either the step did not lead to a decrease in cost or there
- // was numerical failure. In either case we will scale mu up and
- // retry. If it was a numerical failure, we hope that the
- // stronger regularization will make the linear system better
- // conditioned. If it was numerically sane, but there was no
- // decrease in cost, then increasing mu reduces the size of the
- // trust region and we look for a decrease closer to the
- // linearization point.
- ++summary->num_unsuccessful_steps;
- mu = mu * nu;
- nu = 2 * nu;
- }
-
- ++iteration;
-
- total_cost = summary->fixed_cost + cost;
-
- iteration_summary.iteration = iteration;
- iteration_summary.step_is_successful = step_is_successful;
- iteration_summary.cost = total_cost;
- iteration_summary.cost_change = actual_cost_change;
- iteration_summary.gradient_max_norm = gradient_max_norm;
- iteration_summary.step_norm = step_norm;
- iteration_summary.relative_decrease = relative_decrease;
- iteration_summary.mu = mu;
- iteration_summary.eta = options.eta;
- iteration_summary.iteration_time_sec = (time(NULL) - iteration_start_time);
-
- if (options.logging_type >= PER_MINIMIZER_ITERATION) {
- summary->iterations.push_back(iteration_summary);
- }
-
- // Call the various callbacks.
- for (int i = 0; i < options.callbacks.size(); ++i) {
- if (!RunCallback(options.callbacks[i], iteration_summary, summary)) {
- return;
- }
- }
- }
-}
-
-} // namespace internal
-} // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt_strategy.cc b/internal/ceres/levenberg_marquardt_strategy.cc
new file mode 100644
index 0000000..54ff783
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy.cc
@@ -0,0 +1,103 @@
+#include "ceres/levenberg_marquardt_strategy.h"
+
+#include <cmath>
+#include "glog/logging.h"
+#include "ceres/array_utils.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_matrix.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/types.h"
+#include "Eigen/Core"
+
+namespace ceres {
+namespace internal {
+
+LevenbergMarquardtStrategy::LevenbergMarquardtStrategy(
+ const TrustRegionStrategy::Options& options)
+ : linear_solver_(options.linear_solver),
+ radius_(options.initial_radius),
+ max_radius_(options.max_radius),
+ min_diagonal_(options.lm_min_diagonal),
+ max_diagonal_(options.lm_max_diagonal),
+ decrease_factor_(2.0),
+ reuse_diagonal_(false) {
+ CHECK_NOTNULL(linear_solver_);
+ CHECK_GT(min_diagonal_, 0.0);
+ CHECK_LT(min_diagonal_, max_diagonal_);
+ CHECK_GT(max_radius_, 0.0);
+}
+
+LevenbergMarquardtStrategy::~LevenbergMarquardtStrategy() {
+}
+
+LinearSolver::Summary LevenbergMarquardtStrategy::ComputeStep(
+ const TrustRegionStrategy::PerSolveOptions& per_solve_options,
+ SparseMatrix* jacobian,
+ const double* residuals,
+ double* step) {
+ CHECK_NOTNULL(jacobian);
+ CHECK_NOTNULL(residuals);
+ CHECK_NOTNULL(step);
+
+ const int num_parameters = jacobian->num_cols();
+ if (!reuse_diagonal_) {
+ if (diagonal_.rows() != num_parameters) {
+ diagonal_.resize(num_parameters, 1);
+ }
+
+ jacobian->SquaredColumnNorm(diagonal_.data());
+ for (int i = 0; i < num_parameters; ++i) {
+ diagonal_[i] = min(max(diagonal_[i], min_diagonal_), max_diagonal_);
+ }
+ }
+
+ lm_diagonal_ = (diagonal_ / radius_).array().sqrt();
+
+ LinearSolver::PerSolveOptions solve_options;
+ solve_options.D = lm_diagonal_.data();
+ solve_options.q_tolerance = per_solve_options.eta;
+ // Disable r_tolerance checking. Since we only care about
+ // termination via the q_tolerance. As Nash and Sofer show,
+ // r_tolerance based termination is essentially useless in
+ // Truncated Newton methods.
+ solve_options.r_tolerance = -1.0;
+
+ // Invalidate the output array lm_step, so that we can detect if
+ // the linear solver generated numerical garbage. This is known
+ // to happen for the DENSE_QR and then DENSE_SCHUR solver when
+ // the Jacobin is severly rank deficient and mu is too small.
+ InvalidateArray(num_parameters, step);
+ LinearSolver::Summary linear_solver_summary =
+ linear_solver_->Solve(jacobian, residuals, solve_options, step);
+ if (linear_solver_summary.termination_type == FAILURE ||
+ !IsArrayValid(num_parameters, step)) {
+ LOG(WARNING) << "Linear solver failure. Failed to compute a finite step.";
+ linear_solver_summary.termination_type = FAILURE;
+ }
+
+ reuse_diagonal_ = true;
+ return linear_solver_summary;
+}
+
+void LevenbergMarquardtStrategy::StepAccepted(double step_quality) {
+ CHECK_GT(step_quality, 0.0);
+ radius_ = radius_ / std::max(1.0 / 3.0,
+ 1.0 - pow(2.0 * step_quality - 1.0, 3));
+ radius_ = std::min(max_radius_, radius_);
+ decrease_factor_ = 2.0;
+ reuse_diagonal_ = false;
+}
+
+void LevenbergMarquardtStrategy::StepRejected(double step_quality) {
+ radius_ = radius_ / decrease_factor_;
+ decrease_factor_ *= 2.0;
+ reuse_diagonal_ = true;
+}
+
+double LevenbergMarquardtStrategy::Radius() const {
+ return radius_;
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt_strategy.h b/internal/ceres/levenberg_marquardt_strategy.h
new file mode 100644
index 0000000..274e972
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy.h
@@ -0,0 +1,78 @@
+// 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_LEVENBERG_MARQUARDT_STRATEGY_H_
+#define CERES_INTERNAL_LEVENBERG_MARQUARDT_STRATEGY_H_
+
+#include "ceres/linear_solver.h"
+#include "ceres/trust_region_strategy.h"
+
+namespace ceres {
+namespace internal {
+
+// Levenberg-Marquardt step computation and trust region sizing
+// strategy based on on "Methods for Nonlinear Least Squares" by
+// K. Madsen, H.B. Nielsen and O. Tingleff. Available to download from
+//
+// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
+class LevenbergMarquardtStrategy : public TrustRegionStrategy {
+public:
+ LevenbergMarquardtStrategy(const TrustRegionStrategy::Options& options);
+ virtual ~LevenbergMarquardtStrategy();
+
+ // TrustRegionStrategy interface
+ virtual LinearSolver::Summary ComputeStep(
+ const TrustRegionStrategy::PerSolveOptions& per_solve_options,
+ SparseMatrix* jacobian,
+ const double* residuals,
+ double* step);
+ virtual void StepAccepted(double step_quality);
+ virtual void StepRejected(double step_quality);
+ virtual double Radius() const;
+
+ private:
+ LinearSolver* linear_solver_;
+ double radius_;
+ double max_radius_;
+ const double min_diagonal_;
+ const double max_diagonal_;
+ double decrease_factor_;
+ bool reuse_diagonal_;
+ Vector diagonal_; // diagonal_ = diag(J'J)
+ // Scaled copy of diagonal_. Stored here as optimization to prevent
+ // allocations in every iteration and reuse when a step fails and
+ // ComputeStep is called again.
+ Vector lm_diagonal_; // lm_diagonal_ = diagonal_ / radius_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_LEVENBERG_MARQUARDT_STRATEGY_H_
diff --git a/internal/ceres/levenberg_marquardt_strategy_test.cc b/internal/ceres/levenberg_marquardt_strategy_test.cc
new file mode 100644
index 0000000..370d0a5
--- /dev/null
+++ b/internal/ceres/levenberg_marquardt_strategy_test.cc
@@ -0,0 +1,152 @@
+// 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 <glog/logging.h>
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+#include "ceres/mock_log.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/levenberg_marquardt_strategy.h"
+#include "ceres/internal/eigen.h"
+
+using testing::AllOf;
+using testing::AnyNumber;
+using testing::HasSubstr;
+using testing::ScopedMockLog;
+using testing::_;
+
+namespace ceres {
+namespace internal {
+
+const double kTolerance = 1e-16;
+
+// Linear solver that takes as input a vector and checks that the
+// caller passes the same vector as LinearSolver::PerSolveOptions.D.
+class RegularizationCheckingLinearSolver : public DenseSparseMatrixSolver {
+ public:
+ RegularizationCheckingLinearSolver(const int num_cols, const double* diagonal)
+ : num_cols_(num_cols),
+ diagonal_(diagonal) {
+ }
+
+ virtual ~RegularizationCheckingLinearSolver(){}
+
+ private:
+ virtual LinearSolver::Summary SolveImpl(
+ DenseSparseMatrix* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x) {
+ CHECK_NOTNULL(per_solve_options.D);
+ for (int i = 0; i < num_cols_; ++i) {
+ EXPECT_NEAR(per_solve_options.D[i], diagonal_[i], kTolerance)
+ << i << " " << per_solve_options.D[i] << " " << diagonal_[i];
+ }
+ }
+ const int num_cols_;
+ const double* diagonal_;
+};
+
+TEST(LevenbergMarquardtStrategy, AcceptRejectStepRadiusScaling) {
+ TrustRegionStrategy::Options options;
+ options.initial_radius = 2.0;
+ options.max_radius = 20.0;
+ options.lm_min_diagonal = 1e-8;
+ options.lm_max_diagonal = 1e8;
+
+ // We need a non-null pointer here, so anything should do.
+ options.linear_solver = new RegularizationCheckingLinearSolver(0, NULL);
+
+ LevenbergMarquardtStrategy lms(options);
+ EXPECT_EQ(lms.Radius(), options.initial_radius);
+ lms.StepRejected(0.0);
+ EXPECT_EQ(lms.Radius(), 1.0);
+ lms.StepRejected(-1.0);
+ EXPECT_EQ(lms.Radius(), 0.25);
+ lms.StepAccepted(1.0);
+ EXPECT_EQ(lms.Radius(), 0.25 * 3.0);
+ lms.StepAccepted(1.0);
+ EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0);
+ lms.StepAccepted(0.25);
+ EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125);
+ lms.StepAccepted(1.0);
+ EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125 * 3.0);
+ lms.StepAccepted(1.0);
+ EXPECT_EQ(lms.Radius(), 0.25 * 3.0 * 3.0 / 1.125 * 3.0 * 3.0);
+ lms.StepAccepted(1.0);
+ EXPECT_EQ(lms.Radius(), options.max_radius);
+}
+
+TEST(LevenbergMarquardtStrategy, CorrectDiagonalToLinearSolver) {
+ Matrix jacobian(2,3);
+ jacobian.setZero();
+ jacobian(0,0) = 0.0;
+ jacobian(0,1) = 1.0;
+ jacobian(1,1) = 1.0;
+ jacobian(0,2) = 100.0;
+
+ double residual = 1.0;
+ double x[2];
+ DenseSparseMatrix dsm(jacobian);
+
+ TrustRegionStrategy::Options options;
+ options.initial_radius = 2.0;
+ options.max_radius = 20.0;
+ options.lm_min_diagonal = 1e-2;
+ options.lm_max_diagonal = 1e2;
+
+ double diagonal[3];
+ diagonal[0] = options.lm_min_diagonal;
+ diagonal[1] = 2.0;
+ diagonal[2] = options.lm_max_diagonal;
+ for (int i = 0; i < 3; ++i) {
+ diagonal[i] = sqrt(diagonal[i] / options.initial_radius);
+ }
+
+ RegularizationCheckingLinearSolver linear_solver(3, diagonal);
+ options.linear_solver = &linear_solver;
+
+ LevenbergMarquardtStrategy lms(options);
+ TrustRegionStrategy::PerSolveOptions pso;
+
+ {
+ ScopedMockLog log;
+ EXPECT_CALL(log, Log(_, _, _)).Times(AnyNumber());
+ EXPECT_CALL(log, Log(WARNING, _,
+ HasSubstr("Failed to compute a finite step.")));
+
+ LinearSolver::Summary summary = lms.ComputeStep(pso, &dsm, &residual, x);
+ EXPECT_EQ(summary.termination_type, FAILURE);
+ }
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h
index 19b0344..eeda298 100644
--- a/internal/ceres/minimizer.h
+++ b/internal/ceres/minimizer.h
@@ -40,6 +40,8 @@
class Evaluator;
class LinearSolver;
+class SparseMatrix;
+class TrustRegionStrategy;
// Interface for non-linear least squares solvers.
class Minimizer {
@@ -48,55 +50,87 @@
// see solver.h for detailed information about the meaning and
// default values of each of these parameters.
struct Options {
+ Options() {
+ Init(Solver::Options());
+ }
+
explicit Options(const Solver::Options& options) {
+ Init(options);
+ }
+
+ void Init(const Solver::Options& options) {
max_num_iterations = options.max_num_iterations;
max_solver_time_sec = options.max_solver_time_sec;
+ max_step_solver_retries = 5;
gradient_tolerance = options.gradient_tolerance;
parameter_tolerance = options.parameter_tolerance;
function_tolerance = options.function_tolerance;
min_relative_decrease = options.min_relative_decrease;
eta = options.eta;
- tau = options.tau;
- min_mu = options.min_mu;
jacobi_scaling = options.jacobi_scaling;
- crash_and_dump_lsqp_on_failure = options.crash_and_dump_lsqp_on_failure;
lsqp_dump_directory = options.lsqp_dump_directory;
lsqp_iterations_to_dump = options.lsqp_iterations_to_dump;
lsqp_dump_format_type = options.lsqp_dump_format_type;
num_eliminate_blocks = options.num_eliminate_blocks;
- logging_type = options.logging_type;
+ max_num_consecutive_invalid_steps =
+ options.max_num_consecutive_invalid_steps;
+ min_trust_region_radius = options.min_trust_region_radius;
+ evaluator = NULL;
+ trust_region_strategy = NULL;
+ jacobian = NULL;
}
int max_num_iterations;
int max_solver_time_sec;
+
+ // Number of times the linear solver should be retried in case of
+ // numerical failure. The retries are done by exponentially scaling up
+ // mu at each retry. This leads to stronger and stronger
+ // regularization making the linear least squares problem better
+ // conditioned at each retry.
+ int max_step_solver_retries;
double gradient_tolerance;
double parameter_tolerance;
double function_tolerance;
double min_relative_decrease;
double eta;
- double tau;
- double min_mu;
bool jacobi_scaling;
- bool crash_and_dump_lsqp_on_failure;
vector<int> lsqp_iterations_to_dump;
DumpFormatType lsqp_dump_format_type;
string lsqp_dump_directory;
int num_eliminate_blocks;
- LoggingType logging_type;
+ int max_num_consecutive_invalid_steps;
+ int min_trust_region_radius;
// List of callbacks that are executed by the Minimizer at the end
// of each iteration.
//
- // Client owns these pointers.
+ // The Options struct does not own these pointers.
vector<IterationCallback*> callbacks;
+
+ // Object responsible for evaluating the cost, residuals and
+ // Jacobian matrix. The Options struct does not own this pointer.
+ Evaluator* evaluator;
+
+ // Object responsible for actually computing the trust region
+ // step, and sizing the trust region radius. The Options struct
+ // does not own this pointer.
+ TrustRegionStrategy* trust_region_strategy;
+
+ // Object holding the Jacobian matrix. It is assumed that the
+ // sparsity structure of the matrix has already been initialized
+ // and will remain constant for the life time of the
+ // optimization. The Options struct does not own this pointer.
+ SparseMatrix* jacobian;
};
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.
virtual void Minimize(const Options& options,
- Evaluator* evaluator,
- LinearSolver* linear_solver,
- const double* initial_parameters,
- double* final_parameters,
+ double* parameters,
Solver::Summary* summary) = 0;
};
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index d642a54..f5f5d91 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -32,7 +32,6 @@
#include "ceres/solver.h"
#include <vector>
-#include "ceres/levenberg_marquardt.h"
#include "ceres/program.h"
#include "ceres/solver_impl.h"
#include "ceres/stringprintf.h"
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 6f8fa17..4eb3d1f 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -34,18 +34,19 @@
#include <numeric>
#include "ceres/evaluator.h"
#include "ceres/gradient_checking_cost_function.h"
-#include "ceres/levenberg_marquardt.h"
+#include "ceres/iteration_callback.h"
+#include "ceres/levenberg_marquardt_strategy.h"
#include "ceres/linear_solver.h"
#include "ceres/map_util.h"
#include "ceres/minimizer.h"
#include "ceres/parameter_block.h"
+#include "ceres/problem.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/residual_block.h"
#include "ceres/schur_ordering.h"
#include "ceres/stringprintf.h"
-#include "ceres/iteration_callback.h"
-#include "ceres/problem.h"
+#include "ceres/trust_region_minimizer.h"
namespace ceres {
namespace internal {
@@ -104,7 +105,7 @@
summary.gradient_max_norm,
summary.step_norm,
summary.relative_decrease,
- summary.mu,
+ summary.trust_region_radius,
summary.linear_solver_iterations);
if (log_to_stdout_) {
cout << output << endl;
@@ -124,32 +125,40 @@
Program* program,
Evaluator* evaluator,
LinearSolver* linear_solver,
- double* initial_parameters,
- double* final_parameters,
+ double* parameters,
Solver::Summary* summary) {
Minimizer::Options minimizer_options(options);
-
LoggingCallback logging_callback(options.minimizer_progress_to_stdout);
if (options.logging_type != SILENT) {
minimizer_options.callbacks.push_back(&logging_callback);
}
- StateUpdatingCallback updating_callback(program, final_parameters);
+ StateUpdatingCallback updating_callback(program, parameters);
if (options.update_state_every_iteration) {
minimizer_options.callbacks.push_back(&updating_callback);
}
- LevenbergMarquardt levenberg_marquardt;
+ minimizer_options.evaluator = evaluator;
+ scoped_ptr<SparseMatrix> jacobian(evaluator->CreateJacobian());
+ minimizer_options.jacobian = jacobian.get();
- time_t start_minimizer_time_seconds = time(NULL);
- levenberg_marquardt.Minimize(minimizer_options,
- evaluator,
- linear_solver,
- initial_parameters,
- final_parameters,
- summary);
- summary->minimizer_time_in_seconds =
- time(NULL) - start_minimizer_time_seconds;
+ TrustRegionStrategy::Options trust_region_strategy_options;
+ trust_region_strategy_options.linear_solver = linear_solver;
+ trust_region_strategy_options.initial_radius =
+ options.initial_trust_region_radius;
+ trust_region_strategy_options.max_radius = options.max_trust_region_radius;
+ trust_region_strategy_options.lm_min_diagonal = options.lm_min_diagonal;
+ trust_region_strategy_options.lm_max_diagonal = options.lm_max_diagonal;
+ trust_region_strategy_options.trust_region_strategy_type =
+ options.trust_region_strategy_type;
+ scoped_ptr<TrustRegionStrategy> strategy(
+ TrustRegionStrategy::Create(trust_region_strategy_options));
+ minimizer_options.trust_region_strategy = strategy.get();
+
+ TrustRegionMinimizer minimizer;
+ time_t minimizer_start_time = time(NULL);
+ minimizer.Minimize(minimizer_options, parameters, summary);
+ summary->minimizer_time_in_seconds = time(NULL) - minimizer_start_time;
}
void SolverImpl::Solve(const Solver::Options& original_options,
@@ -259,19 +268,17 @@
}
// The optimizer works on contiguous parameter vectors; allocate some.
- Vector initial_parameters(reduced_program->NumParameters());
- Vector optimized_parameters(reduced_program->NumParameters());
+ Vector parameters(reduced_program->NumParameters());
// Collect the discontiguous parameters into a contiguous state vector.
- reduced_program->ParameterBlocksToStateVector(&initial_parameters[0]);
+ reduced_program->ParameterBlocksToStateVector(parameters.data());
// Run the optimization.
Minimize(options,
reduced_program.get(),
evaluator.get(),
linear_solver.get(),
- initial_parameters.data(),
- optimized_parameters.data(),
+ parameters.data(),
summary);
// If the user aborted mid-optimization or the optimization
@@ -283,7 +290,7 @@
}
// Push the contiguous optimized parameters back to the user's parameters.
- reduced_program->StateVectorToParameterBlocks(&optimized_parameters[0]);
+ reduced_program->StateVectorToParameterBlocks(parameters.data());
reduced_program->CopyParameterBlockStateToUserState();
// Return the final cost and residuals for the original problem.
diff --git a/internal/ceres/solver_impl.h b/internal/ceres/solver_impl.h
index 957ebcc..7dee03c 100644
--- a/internal/ceres/solver_impl.h
+++ b/internal/ceres/solver_impl.h
@@ -92,8 +92,7 @@
Program* program,
Evaluator* evaluator,
LinearSolver* linear_solver,
- double* initial_parameters,
- double* final_parameters,
+ double* parameters,
Solver::Summary* summary);
// Remove the fixed or unused parameter blocks and residuals
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
new file mode 100644
index 0000000..bb6c4a1
--- /dev/null
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -0,0 +1,445 @@
+// 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/trust_region_minimizer.h"
+
+#include <algorithm>
+#include <cstdlib>
+#include <cmath>
+#include <cstring>
+#include <string>
+#include <vector>
+#include <glog/logging.h>
+
+#include "Eigen/Core"
+#include "ceres/array_utils.h"
+#include "ceres/evaluator.h"
+#include "ceres/linear_least_squares_problems.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/sparse_matrix.h"
+#include "ceres/trust_region_strategy.h"
+#include "ceres/types.h"
+
+namespace ceres {
+namespace internal {
+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 TrustRegionMinimizer::RunCallbacks(
+ const Minimizer::Options& options_,
+ 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,
+ double* scale) const {
+ jacobian.SquaredColumnNorm(scale);
+ for (int i = 0; i < jacobian.num_cols(); ++i) {
+ scale[i] = 1.0 / (kEpsilon + sqrt(scale[i]));
+ }
+}
+
+void TrustRegionMinimizer::Init(const Minimizer::Options& options) {
+ options_ = options;
+ sort(options_.lsqp_iterations_to_dump.begin(),
+ options_.lsqp_iterations_to_dump.end());
+}
+
+bool TrustRegionMinimizer::MaybeDumpLinearLeastSquaresProblem(
+ const int iteration,
+ const SparseMatrix* jacobian,
+ const double* residuals,
+ const double* step) const {
+ // TODO(sameeragarwal): Since the use of trust_region_radius has
+ // moved inside TrustRegionStrategy, its not clear how we dump the
+ // regularization vector/matrix anymore.
+ //
+ // Doing this right requires either an API change to the
+ // TrustRegionStrategy and/or how LinearLeastSquares problems are
+ // stored on disk.
+ //
+ // For now, we will just not dump the regularizer.
+ return (!binary_search(options_.lsqp_iterations_to_dump.begin(),
+ options_.lsqp_iterations_to_dump.end(),
+ iteration) ||
+ DumpLinearLeastSquaresProblem(options_.lsqp_dump_directory,
+ iteration,
+ options_.lsqp_dump_format_type,
+ jacobian,
+ NULL,
+ residuals,
+ step,
+ options_.num_eliminate_blocks));
+}
+
+void TrustRegionMinimizer::Minimize(const Minimizer::Options& options,
+ double* parameters,
+ Solver::Summary* summary) {
+ time_t start_time = time(NULL);
+ time_t iteration_start_time = start_time;
+ Init(options);
+
+ summary->termination_type = NO_CONVERGENCE;
+ summary->num_successful_steps = 0;
+ summary->num_unsuccessful_steps = 0;
+
+ Evaluator* evaluator = CHECK_NOTNULL(options_.evaluator);
+ SparseMatrix* jacobian = CHECK_NOTNULL(options_.jacobian);
+ TrustRegionStrategy* strategy = CHECK_NOTNULL(options_.trust_region_strategy);
+
+ const int num_parameters = evaluator->NumParameters();
+ const int num_effective_parameters = evaluator->NumEffectiveParameters();
+ const int num_residuals = evaluator->NumResiduals();
+
+ VectorRef x(parameters, num_parameters);
+ double x_norm = x.norm();
+
+ Vector residuals(num_residuals);
+ Vector trust_region_step(num_effective_parameters);
+ Vector delta(num_effective_parameters);
+ Vector x_plus_delta(num_parameters);
+ Vector gradient(num_effective_parameters);
+ Vector model_residuals(num_residuals);
+ Vector scale(num_effective_parameters);
+
+ IterationSummary iteration_summary;
+ iteration_summary.iteration = 0;
+ iteration_summary.step_is_valid=false;
+ iteration_summary.step_is_successful=false;
+ iteration_summary.cost = summary->initial_cost;
+ iteration_summary.cost_change = 0.0;
+ iteration_summary.gradient_max_norm = 0.0;
+ iteration_summary.step_norm = 0.0;
+ iteration_summary.relative_decrease = 0.0;
+ iteration_summary.trust_region_radius = strategy->Radius();
+ // TODO(sameeragarwal): Rename eta to linear_solver_accuracy or
+ // something similar across the board.
+ iteration_summary.eta = options_.eta;
+ iteration_summary.linear_solver_iterations = 0;
+ iteration_summary.step_solver_time_sec = 0;
+
+ // Do initial cost and Jacobian evaluation.
+ double cost = 0.0;
+ if (!evaluator->Evaluate(x.data(), &cost, residuals.data(), jacobian)) {
+ LOG(WARNING) << "Terminating: Residual and Jacobian evaluation failed.";
+ summary->termination_type = NUMERICAL_FAILURE;
+ return;
+ }
+
+ // Compute the fixed part of the cost.
+ //
+ // This is a poor way to do this computation. Even if fixed_cost is
+ // zero, because we are subtracting two possibly large numbers, we
+ // are depending on exact cancellation to give us a zero here. But
+ // initial_cost and cost have been computed by two different
+ // evaluators. One which runs on the whole problem (in
+ // solver_impl.cc) in single threaded mode and another which runs
+ // here on the reduced problem, so fixed_cost can (and does) contain
+ // some numerical garbage with a relative magnitude of 1e-14.
+ //
+ // The right way to do this, would be to compute the fixed cost on
+ // just the set of residual blocks which are held constant and were
+ // removed from the original problem when the reduced problem was
+ // constructed.
+ summary->fixed_cost = summary->initial_cost - cost;
+
+ gradient.setZero();
+ jacobian->LeftMultiply(residuals.data(), gradient.data());
+ iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
+
+ if (options_.jacobi_scaling) {
+ EstimateScale(*jacobian, scale.data());
+ jacobian->ScaleColumns(scale.data());
+ } else {
+ scale.setOnes();
+ }
+
+ const double absolute_gradient_tolerance =
+ options_.gradient_tolerance *
+ max(iteration_summary.gradient_max_norm, kEpsilon);
+ 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 /
+ summary->iterations[0].gradient_max_norm
+ << " <= " << options_.gradient_tolerance;
+ return;
+ }
+
+ // Call the various callbacks.
+ iteration_summary.iteration_time_sec = time(NULL) - iteration_start_time;
+ summary->iterations.push_back(iteration_summary);
+ switch (RunCallbacks(options_, 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;
+ while (true) {
+ iteration_start_time = time(NULL);
+ if (iteration_summary.iteration >= options_.max_num_iterations) {
+ summary->termination_type = NO_CONVERGENCE;
+ VLOG(1) << "Terminating: Maximum number of iterations reached.";
+ break;
+ }
+
+ if ((start_time - iteration_start_time) >= options_.max_solver_time_sec) {
+ summary->termination_type = NO_CONVERGENCE;
+ VLOG(1) << "Terminating: Maximum solver time reached.";
+ break;
+ }
+
+ iteration_summary = IterationSummary();
+ iteration_summary = summary->iterations.back();
+ iteration_summary.iteration = summary->iterations.back().iteration + 1;
+ iteration_summary.step_is_valid = false;
+ iteration_summary.step_is_successful = false;
+
+ const time_t strategy_start_time = time(NULL);
+ TrustRegionStrategy::PerSolveOptions per_solve_options;
+ per_solve_options.eta = options_.eta;
+ LinearSolver::Summary strategy_summary =
+ strategy->ComputeStep(per_solve_options,
+ jacobian,
+ residuals.data(),
+ trust_region_step.data());
+
+ iteration_summary.step_solver_time_sec = time(NULL) - strategy_start_time;
+ iteration_summary.linear_solver_iterations =
+ strategy_summary.num_iterations;
+
+ if (!MaybeDumpLinearLeastSquaresProblem(iteration_summary.iteration,
+ jacobian,
+ residuals.data(),
+ trust_region_step.data())) {
+ LOG(FATAL) << "Tried writing linear least squares problem: "
+ << options.lsqp_dump_directory << "but failed.";
+ }
+
+ double new_model_cost = 0.0;
+ if (strategy_summary.termination_type != FAILURE) {
+ // new_model_cost = 1/2 |J * step - f|^2
+ model_residuals = -residuals;
+ jacobian->RightMultiply(trust_region_step.data(), model_residuals.data());
+ new_model_cost = model_residuals.squaredNorm() / 2.0;
+
+ // In exact arithmetic, this would never be the case. But poorly
+ // conditioned matrices can give rise to situations where the
+ // new_model_cost can actually be larger than half the squared
+ // norm of the residual vector. We allow for small tolerance
+ // around cost and beyond that declare the step to be invalid.
+ if (cost < (new_model_cost - kEpsilon)) {
+ VLOG(1) << "Invalid step: current_cost: " << cost
+ << " new_model_cost " << new_model_cost;
+ } else {
+ iteration_summary.step_is_valid = true;
+ }
+ }
+
+ if (!iteration_summary.step_is_valid) {
+ // Invalid steps can happen due to a number of reasons, and we
+ // allow a limited number of successive failures, and return with
+ // NUMERICAL_FAILURE if this limit is exceeded.
+ if (++num_consecutive_invalid_steps >=
+ options_.max_num_consecutive_invalid_steps) {
+ summary->termination_type = NUMERICAL_FAILURE;
+ LOG(WARNING) << "Terminating. Number of successive invalid steps more "
+ << "than "
+ << "Solver::Options::max_num_consecutive_invalid_steps: "
+ << options_.max_num_consecutive_invalid_steps;
+ return;
+ }
+
+ // We are going to try and reduce the trust region radius and
+ // solve again. To do this, we are going to treat this iteration
+ // as an unsuccessful iteration. Since the various callbacks are
+ // still executed, we are going to fill the iteration summary
+ // with data that assumes a step of length zero and no progress.
+ iteration_summary.cost = cost;
+ iteration_summary.cost_change = 0.0;
+ iteration_summary.gradient_max_norm =
+ summary->iterations.back().gradient_max_norm;
+ iteration_summary.step_norm = 0.0;
+ iteration_summary.relative_decrease = 0.0;
+ iteration_summary.eta = options_.eta;
+ } else {
+ // The step is numerically valid, so now we can judge its quality.
+ num_consecutive_invalid_steps = 0;
+
+ // We allow some slop around 0, and clamp the model_cost_change
+ // at kEpsilon from below.
+ //
+ // There is probably a better way to do this, as it is going to
+ // create problems for problems where the objective function is
+ // kEpsilon close to zero.
+ const double model_cost_change = max(kEpsilon, cost - new_model_cost);
+
+ // Undo the Jacobian column scaling.
+ delta = -(trust_region_step.array() * scale.array()).matrix();
+ iteration_summary.step_norm = delta.norm();
+
+ // Convergence based on parameter_tolerance.
+ const double step_size_tolerance = options_.parameter_tolerance *
+ (x_norm + options_.parameter_tolerance);
+ if (iteration_summary.step_norm <= step_size_tolerance) {
+ VLOG(1) << "Terminating. Parameter tolerance reached. "
+ << "relative step_norm: "
+ << iteration_summary.step_norm /
+ (x_norm + options_.parameter_tolerance)
+ << " <= " << options_.parameter_tolerance;
+ summary->termination_type = PARAMETER_TOLERANCE;
+ return;
+ }
+
+ if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
+ summary->termination_type = NUMERICAL_FAILURE;
+ LOG(WARNING) << "Terminating. Failed to compute "
+ << "Plus(x, delta, x_plus_delta).";
+ return;
+ }
+
+ // Try this step.
+ double new_cost;
+ if (!evaluator->Evaluate(x_plus_delta.data(), &new_cost, NULL, NULL)) {
+ summary->termination_type = NUMERICAL_FAILURE;
+ LOG(WARNING) << "Terminating: Cost evaluation failed.";
+ return;
+ }
+
+ VLOG(2) << "old cost: " << cost << "new cost: " << new_cost;
+ iteration_summary.cost_change = cost - new_cost;
+ const double absolute_function_tolerance =
+ options_.function_tolerance * cost;
+ if (fabs(iteration_summary.cost_change) < absolute_function_tolerance) {
+ VLOG(1) << "Terminating. Function tolerance reached. "
+ << "|cost_change|/cost: " << fabs(iteration_summary.cost_change) / cost
+ << " <= " << options_.function_tolerance;
+ summary->termination_type = FUNCTION_TOLERANCE;
+ return;
+ }
+
+ iteration_summary.relative_decrease =
+ iteration_summary.cost_change / model_cost_change;
+ iteration_summary.step_is_successful =
+ iteration_summary.relative_decrease > options_.min_relative_decrease;
+ }
+
+ if (iteration_summary.step_is_successful) {
+ ++summary->num_successful_steps;
+ strategy->StepAccepted(iteration_summary.relative_decrease);
+ x = x_plus_delta;
+ x_norm = x.norm();
+ // Step looks good, evaluate the residuals and Jacobian at this
+ // point.
+ if (!evaluator->Evaluate(x.data(), &cost, residuals.data(), jacobian)) {
+ summary->termination_type = NUMERICAL_FAILURE;
+ LOG(WARNING) << "Terminating: Residual and Jacobian evaluation failed.";
+ return;
+ }
+
+ gradient.setZero();
+ jacobian->LeftMultiply(residuals.data(), gradient.data());
+ iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>();
+
+ 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 /
+ summary->iterations[0].gradient_max_norm
+ << " <= " << options_.gradient_tolerance;
+ return;
+ }
+
+ if (options_.jacobi_scaling) {
+ jacobian->ScaleColumns(scale.data());
+ }
+ } else {
+ ++summary->num_unsuccessful_steps;
+ strategy->StepRejected(iteration_summary.relative_decrease);
+ }
+
+ iteration_summary.cost = cost + summary->fixed_cost;
+ iteration_summary.trust_region_radius = strategy->Radius();
+ if (iteration_summary.trust_region_radius <
+ options_.min_trust_region_radius) {
+ summary->termination_type = PARAMETER_TOLERANCE;
+ VLOG(1) << "Termination. Minimum trust region radius reached.";
+ return;
+ }
+
+ iteration_summary.iteration_time_sec = time(NULL) - iteration_start_time;
+ summary->iterations.push_back(iteration_summary);
+ switch (RunCallbacks(options_, 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";
+ }
+ }
+}
+
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/levenberg_marquardt.h b/internal/ceres/trust_region_minimizer.h
similarity index 62%
rename from internal/ceres/levenberg_marquardt.h
rename to internal/ceres/trust_region_minimizer.h
index d00bb90..4337b18 100644
--- a/internal/ceres/levenberg_marquardt.h
+++ b/internal/ceres/trust_region_minimizer.h
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
+// 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
@@ -27,39 +27,42 @@
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sameeragarwal@google.com (Sameer Agarwal)
-//
-// Implmentation of Levenberg Marquardt algorithm based on "Methods for
-// Nonlinear Least Squares" by K. Madsen, H.B. Nielsen and
-// O. Tingleff. Available to download from
-//
-// http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
-//
-#ifndef CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
-#define CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
+#ifndef CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
+#define CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
#include "ceres/minimizer.h"
#include "ceres/solver.h"
+#include "ceres/types.h"
namespace ceres {
namespace internal {
-class Evaluator;
-class LinearSolver;
-
-class LevenbergMarquardt : public Minimizer {
+// Generic trust region minimization algorithm. The heavy lifting is
+// done by a TrustRegionStrategy object passed in as one of the
+// arguments to the Minimize method.
+//
+// For example usage, see SolverImpl::Minimize.
+class TrustRegionMinimizer : public Minimizer {
public:
- virtual ~LevenbergMarquardt();
-
+ ~TrustRegionMinimizer() {}
virtual void Minimize(const Minimizer::Options& options,
- Evaluator* evaluator,
- LinearSolver* linear_solver,
- const double* initial_parameters,
- double* final_parameters,
+ double* parameters,
Solver::Summary* summary);
+
+ private:
+ void Init(const Minimizer::Options& options);
+ void EstimateScale(const SparseMatrix& jacobian, double* scale) const;
+ CallbackReturnType RunCallbacks(const Minimizer::Options& options,
+ const IterationSummary& iteration_summary);
+ bool MaybeDumpLinearLeastSquaresProblem( const int iteration,
+ const SparseMatrix* jacobian,
+ const double* residuals,
+ const double* step) const;
+
+ Minimizer::Options options_;
};
} // namespace internal
} // namespace ceres
-
-#endif // CERES_INTERNAL_LEVENBERG_MARQUARDT_H_
+#endif // CERES_INTERNAL_TRUST_REGION_MINIMIZER_H_
diff --git a/internal/ceres/levenberg_marquardt_test.cc b/internal/ceres/trust_region_minimizer_test.cc
similarity index 72%
rename from internal/ceres/levenberg_marquardt_test.cc
rename to internal/ceres/trust_region_minimizer_test.cc
index 020abfa..b4ef601 100644
--- a/internal/ceres/levenberg_marquardt_test.cc
+++ b/internal/ceres/trust_region_minimizer_test.cc
@@ -1,5 +1,5 @@
// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
+// 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
@@ -27,19 +27,21 @@
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
+// sameeragarwal@google.com (Sameer Agarwal)
//
-// This tests the Levenberg-Marquardt loop using a direct Evaluator
-// implementation, rather than having a test that goes through all the Program
-// and Problem machinery.
+// This tests the TrustRegionMinimizer loop using a direct Evaluator
+// implementation, rather than having a test that goes through all the
+// Program and Problem machinery.
#include <cmath>
#include "ceres/dense_qr_solver.h"
#include "ceres/dense_sparse_matrix.h"
#include "ceres/evaluator.h"
-#include "ceres/levenberg_marquardt.h"
+#include "ceres/internal/port.h"
+#include "ceres/levenberg_marquardt_strategy.h"
#include "ceres/linear_solver.h"
#include "ceres/minimizer.h"
-#include "ceres/internal/port.h"
+#include "ceres/trust_region_minimizer.h"
#include "gtest/gtest.h"
namespace ceres {
@@ -181,67 +183,74 @@
// Templated function to hold a subset of the columns fixed and check
// if the solver converges to the optimal values or not.
template<bool col1, bool col2, bool col3, bool col4>
-void IsSolveSuccessful() {
- LevenbergMarquardt lm;
+void IsLevenbergMarquardtSolveSuccessful() {
Solver::Options solver_options;
+ LinearSolver::Options linear_solver_options;
+ DenseQRSolver linear_solver(linear_solver_options);
+
+ double parameters[4] = { 3, -1, 0, 1.0 };
+
+ // If the column is inactive, then set its value to the optimal
+ // value.
+ parameters[0] = (col1 ? parameters[0] : 0.0);
+ parameters[1] = (col2 ? parameters[1] : 0.0);
+ parameters[2] = (col3 ? parameters[2] : 0.0);
+ parameters[3] = (col4 ? parameters[3] : 0.0);
+
+ PowellEvaluator2<col1, col2, col3, col4> powell_evaluator;
+ scoped_ptr<SparseMatrix> jacobian(powell_evaluator.CreateJacobian());
+
Minimizer::Options minimizer_options(solver_options);
minimizer_options.gradient_tolerance = 1e-26;
minimizer_options.function_tolerance = 1e-26;
minimizer_options.parameter_tolerance = 1e-26;
- LinearSolver::Options linear_solver_options;
- DenseQRSolver linear_solver(linear_solver_options);
+ minimizer_options.evaluator = &powell_evaluator;
+ minimizer_options.jacobian = jacobian.get();
- double initial_parameters[4] = { 3, -1, 0, 1.0 };
- double final_parameters[4] = { -1.0, -1.0, -1.0, -1.0 };
+ TrustRegionStrategy::Options trust_region_strategy_options;
+ trust_region_strategy_options.linear_solver = &linear_solver;
+ trust_region_strategy_options.initial_radius = 1e4;
+ trust_region_strategy_options.max_radius = 1e20;
+ trust_region_strategy_options.lm_min_diagonal = 1e-6;
+ trust_region_strategy_options.lm_max_diagonal = 1e32;
+ scoped_ptr<TrustRegionStrategy> strategy(
+ TrustRegionStrategy::Create(trust_region_strategy_options));
+ minimizer_options.trust_region_strategy = strategy.get();
- // If the column is inactive, then set its value to the optimal
- // value.
- initial_parameters[0] = (col1 ? initial_parameters[0] : 0.0);
- initial_parameters[1] = (col2 ? initial_parameters[1] : 0.0);
- initial_parameters[2] = (col3 ? initial_parameters[2] : 0.0);
- initial_parameters[3] = (col4 ? initial_parameters[3] : 0.0);
-
- PowellEvaluator2<col1, col2, col3, col4> powell_evaluator;
-
+ TrustRegionMinimizer minimizer;
Solver::Summary summary;
- lm.Minimize(minimizer_options,
- &powell_evaluator,
- &linear_solver,
- initial_parameters,
- final_parameters,
- &summary);
+ minimizer.Minimize(minimizer_options, parameters, &summary);
// The minimum is at x1 = x2 = x3 = x4 = 0.
- EXPECT_NEAR(0.0, final_parameters[0], 0.001);
- EXPECT_NEAR(0.0, final_parameters[1], 0.001);
- EXPECT_NEAR(0.0, final_parameters[2], 0.001);
- EXPECT_NEAR(0.0, final_parameters[3], 0.001);
+ EXPECT_NEAR(0.0, parameters[0], 0.001);
+ EXPECT_NEAR(0.0, parameters[1], 0.001);
+ EXPECT_NEAR(0.0, parameters[2], 0.001);
+ EXPECT_NEAR(0.0, parameters[3], 0.001);
};
-TEST(LevenbergMarquardt, PowellsSingularFunction) {
+TEST(TrustRegionMinimizer, PowellsSingularFunction) {
// This case is excluded because this has a local minimum and does
// not find the optimum. This should not affect the correctness of
// this test since we are testing all the other 14 combinations of
// column activations.
+ //
+ // IsSolveSuccessful<true, true, false, true>();
- // IsSolveSuccessful<true, true, false, true>();
-
- IsSolveSuccessful<true, true, true, true>();
- IsSolveSuccessful<true, true, true, false>();
- IsSolveSuccessful<true, false, true, true>();
- IsSolveSuccessful<false, true, true, true>();
- IsSolveSuccessful<true, true, false, false>();
- IsSolveSuccessful<true, false, true, false>();
- IsSolveSuccessful<false, true, true, false>();
- IsSolveSuccessful<true, false, false, true>();
- IsSolveSuccessful<false, true, false, true>();
- IsSolveSuccessful<false, false, true, true>();
- IsSolveSuccessful<true, false, false, false>();
- IsSolveSuccessful<false, true, false, false>();
- IsSolveSuccessful<false, false, true, false>();
- IsSolveSuccessful<false, false, false, true>();
+ IsLevenbergMarquardtSolveSuccessful<true, true, true, true>();
+ IsLevenbergMarquardtSolveSuccessful<true, true, true, false>();
+ IsLevenbergMarquardtSolveSuccessful<true, false, true, true>();
+ IsLevenbergMarquardtSolveSuccessful<false, true, true, true>();
+ IsLevenbergMarquardtSolveSuccessful<true, true, false, false>();
+ IsLevenbergMarquardtSolveSuccessful<true, false, true, false>();
+ IsLevenbergMarquardtSolveSuccessful<false, true, true, false>();
+ IsLevenbergMarquardtSolveSuccessful<true, false, false, true>();
+ IsLevenbergMarquardtSolveSuccessful<false, true, false, true>();
+ IsLevenbergMarquardtSolveSuccessful<false, false, true, true>();
+ IsLevenbergMarquardtSolveSuccessful<true, false, false, false>();
+ IsLevenbergMarquardtSolveSuccessful<false, true, false, false>();
+ IsLevenbergMarquardtSolveSuccessful<false, false, true, false>();
+ IsLevenbergMarquardtSolveSuccessful<false, false, false, true>();
}
-
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/trust_region_strategy.cc b/internal/ceres/trust_region_strategy.cc
new file mode 100644
index 0000000..5c69bd4
--- /dev/null
+++ b/internal/ceres/trust_region_strategy.cc
@@ -0,0 +1,23 @@
+#include "ceres/trust_region_strategy.h"
+#include "ceres/levenberg_marquardt_strategy.h"
+
+namespace ceres {
+namespace internal {
+
+TrustRegionStrategy::~TrustRegionStrategy() {}
+
+TrustRegionStrategy* TrustRegionStrategy::Create(const Options& options) {
+ switch (options.trust_region_strategy_type) {
+ case LEVENBERG_MARQUARDT:
+ return new LevenbergMarquardtStrategy(options);
+ default:
+ LOG(FATAL) << "Unknown trust region strategy: "
+ << options.trust_region_strategy_type;
+ }
+
+ LOG(FATAL) << "Unknown trust region strategy: "
+ << options.trust_region_strategy_type;
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/trust_region_strategy.h b/internal/ceres/trust_region_strategy.h
new file mode 100644
index 0000000..afa4400
--- /dev/null
+++ b/internal/ceres/trust_region_strategy.h
@@ -0,0 +1,109 @@
+// 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_TRUST_REGION_STRATEGY_H_
+#define CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
+
+#include "ceres/linear_solver.h"
+
+namespace ceres {
+namespace internal {
+
+// Interface for classes implementing various trust region strategies
+// for nonlinear least squares problems.
+//
+// The object is expected to maintain and update a trust region
+// radius, which it then uses to solve for the trust region step using
+// the jacobian matrix and residual vector.
+//
+// Here the term trust region radius is used loosely, as the strategy
+// is free to treat it as guidance and violate it as need be. e.g.,
+// the LevenbergMarquardtStrategy uses the inverse of the trust region
+// radius to scale the damping term, which controls the step size, but
+// does not set a hard limit on its size.
+class TrustRegionStrategy {
+public:
+ struct Options {
+ Options()
+ : trust_region_strategy_type(LEVENBERG_MARQUARDT),
+ initial_radius(1e4),
+ max_radius(1e32),
+ lm_min_diagonal(1e-6),
+ lm_max_diagonal(1e32) {
+ }
+
+ TrustRegionStrategyType trust_region_strategy_type;
+ // Linear solver used for actually solving the trust region step.
+ LinearSolver* linear_solver;
+ double initial_radius;
+ double max_radius;
+
+ // Minimum and maximum values of the diagonal damping matrix used
+ // by LevenbergMarquardtStrategy.
+ double lm_min_diagonal;
+ double lm_max_diagonal;
+ };
+
+ // Per solve options.
+ struct PerSolveOptions {
+ // Forcing sequence for inexact solves.
+ double eta;
+ };
+
+ ~TrustRegionStrategy();
+
+ // Use the current radius to solve for the trust region step.
+ virtual LinearSolver::Summary ComputeStep(
+ const PerSolveOptions& per_solve_options,
+ SparseMatrix* jacobian,
+ const double* residuals,
+ double* step) = 0;
+
+ // Inform the strategy that the current step has been accepted, and
+ // that the ratio of the decrease in the non-linear objective to the
+ // decrease in the trust region model is step_quality.
+ virtual void StepAccepted(double step_quality) = 0;
+
+ // Inform the strategy that the current step has been rejected, and
+ // that the ratio of the decrease in the non-linear objective to the
+ // decrease in the trust region model is step_quality.
+ virtual void StepRejected(double step_quality) = 0;
+
+ // Current trust region radius.
+ virtual double Radius() const = 0;
+
+ // Factory.
+ static TrustRegionStrategy* Create(const Options& options);
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
