Non-monotonic trust region algorithm. Non-monotonic trust region algorithm based on the work of Phil Toint, as described in Non-monotone trust region algorithms for nonlinear optimization subject to convex constraints. Philippe L. Toint Mathematical Programming 77 (1997), 69-94. Change-Id: I199ecc644e8d1a8cb43666052aef66fb93e15569
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 5e1b3a6..bb8d3b4 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc
@@ -86,6 +86,8 @@ DEFINE_bool(use_block_amd, true, "Use a block oriented fill reducing ordering."); DEFINE_string(trust_region_strategy, "lm", "Options are: lm, dogleg"); DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds."); +DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use" + " nonmonotic steps"); namespace ceres { namespace examples { @@ -221,6 +223,7 @@ options->num_threads = FLAGS_num_threads; options->eta = FLAGS_eta; options->max_solver_time_in_seconds = FLAGS_max_solver_time; + options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps; if (FLAGS_trust_region_strategy == "lm") { options->trust_region_strategy_type = LEVENBERG_MARQUARDT; } else if (FLAGS_trust_region_strategy == "dogleg") {
diff --git a/include/ceres/solver.h b/include/ceres/solver.h index 4d59bde..8a58cb1 100644 --- a/include/ceres/solver.h +++ b/include/ceres/solver.h
@@ -58,6 +58,8 @@ // Default constructor that sets up a generic sparse problem. Options() { trust_region_strategy_type = LEVENBERG_MARQUARDT; + use_nonmonotonic_steps = false; + max_consecutive_nonmonotonic_steps = 5; max_num_iterations = 50; max_solver_time_in_seconds = 1e9; num_threads = 1; @@ -119,6 +121,34 @@ TrustRegionStrategyType trust_region_strategy_type; + // The classical trust region methods are descent methods, in that + // they only accept a point if it strictly reduces the value of + // the objective function. + // + // Relaxing this requirement allows the algorithm to be more + // efficient in the long term at the cost of some local increase + // in the value of the objective function. + // + // This is because allowing for non-decreasing objective function + // values in a princpled manner allows the algorithm to "jump over + // boulders" as the method is not restricted to move into narrow + // valleys while preserving its convergence properties. + // + // Setting use_nonmonotonic_steps to true enables the + // non-monotonic trust region algorithm as described by Conn, + // Gould & Toint in "Trust Region Methods", Section 10.1. + // + // The parameter max_consecutive_nonmonotonic_steps controls the + // window size used by the step selection algorithm to accept + // non-monotonic steps. + // + // Even though the value of the objective function may be larger + // than the minimum value encountered over the course of the + // optimization, the final parameters returned to the user are the + // ones corresponding to the minimum cost over all iterations. + bool use_nonmonotonic_steps; + int max_consecutive_nonmonotonic_steps; + // Maximum number of iterations for the minimizer to run for. int max_num_iterations;
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h index 70b530f..28b76ce 100644 --- a/internal/ceres/minimizer.h +++ b/internal/ceres/minimizer.h
@@ -68,6 +68,9 @@ min_relative_decrease = options.min_relative_decrease; eta = options.eta; jacobi_scaling = options.jacobi_scaling; + use_nonmonotonic_steps = options.use_nonmonotonic_steps; + max_consecutive_nonmonotonic_steps = + options.max_consecutive_nonmonotonic_steps; 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; @@ -96,6 +99,8 @@ double min_relative_decrease; double eta; bool jacobi_scaling; + bool use_nonmonotonic_steps; + bool max_consecutive_nonmonotonic_steps; vector<int> lsqp_iterations_to_dump; DumpFormatType lsqp_dump_format_type; string lsqp_dump_directory;
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc index 357d822..a6cb830 100644 --- a/internal/ceres/trust_region_minimizer.cc +++ b/internal/ceres/trust_region_minimizer.cc
@@ -132,7 +132,8 @@ const int num_effective_parameters = evaluator->NumEffectiveParameters(); const int num_residuals = evaluator->NumResiduals(); - VectorRef x(parameters, num_parameters); + VectorRef x_min(parameters, num_parameters); + Vector x = x_min; double x_norm = x.norm(); Vector residuals(num_residuals); @@ -145,8 +146,8 @@ IterationSummary iteration_summary; iteration_summary.iteration = 0; - iteration_summary.step_is_valid=false; - iteration_summary.step_is_successful=false; + 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; @@ -167,6 +168,13 @@ return; } + int num_consecutive_nonmonotonic_steps = 0; + double minimum_cost = cost; + double reference_cost = cost; + double accumulated_reference_model_cost_change = 0.0; + double candidate_cost = cost; + double accumulated_candidate_model_cost_change = 0.0; + gradient.setZero(); jacobian->LeftMultiply(residuals.data(), gradient.data()); iteration_summary.gradient_max_norm = gradient.lpNorm<Eigen::Infinity>(); @@ -366,10 +374,43 @@ return; } - iteration_summary.relative_decrease = + const double relative_decrease = iteration_summary.cost_change / model_cost_change; + + const double historical_relative_decrease = + (reference_cost - new_cost) / + (accumulated_reference_model_cost_change + model_cost_change); + + // If monotonic steps are being used, then the relative_decrease + // is the usual ratio of the change in objective function value + // divided by the change in model cost. + // + // If non-monotonic steps are allowed, then we take the maximum + // of the relative_decrease and the + // historical_relative_decrease, which measures the increase + // from a reference iteration. The model cost change is + // estimated by accumulating the model cost changes since the + // reference iteration. The historical relative_decrease offers + // a boost to a step which is not too bad compared to the + // reference iteration, allowing for non-monotonic steps. + iteration_summary.relative_decrease = + options.use_nonmonotonic_steps + ? max(relative_decrease, historical_relative_decrease) + : relative_decrease; + iteration_summary.step_is_successful = iteration_summary.relative_decrease > options_.min_relative_decrease; + + if (iteration_summary.step_is_successful) { + accumulated_candidate_model_cost_change += model_cost_change; + accumulated_reference_model_cost_change += model_cost_change; + if (relative_decrease <= options_.min_relative_decrease) { + VLOG(2) << "Non-monotonic step! " + << " relative_decrease: " << relative_decrease + << " historical_relative_decrease: " + << historical_relative_decrease; + } + } } if (iteration_summary.step_is_successful) { @@ -377,6 +418,7 @@ 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(), @@ -405,6 +447,47 @@ if (options_.jacobi_scaling) { jacobian->ScaleColumns(scale.data()); } + + // Update the best, reference and candidate iterates. + // + // Based on algorithm 10.1.2 (page 357) of "Trust Region + // Methods" by Conn Gould & Toint, or equations 33-40 of + // "Non-monotone trust-region algorithms for nonlinear + // optimization subject to convex constraints" by Phil Toint, + // Mathematical Programming, 77, 1997. + if (cost < minimum_cost) { + // A step that improves solution quality was found. + x_min = x; + minimum_cost = cost; + // Set the candidate iterate to the current point. + candidate_cost = cost; + num_consecutive_nonmonotonic_steps = 0; + accumulated_candidate_model_cost_change = 0.0; + } else { + ++num_consecutive_nonmonotonic_steps; + if (cost > candidate_cost) { + // The current iterate is has a higher cost than the + // candidate iterate. Set the candidate to this point. + VLOG(2) << "Updating the candidate iterate to the current point."; + candidate_cost = cost; + accumulated_candidate_model_cost_change = 0.0; + } + + // At this point we have made too many non-monotonic steps and + // we are going to reset the value of the reference iterate so + // as to force the algorithm to descend. + // + // This is the case because the candidate iterate has a value + // greater than minimum_cost but smaller than the reference + // iterate. + if (num_consecutive_nonmonotonic_steps == + options.max_consecutive_nonmonotonic_steps) { + VLOG(2) << "Resetting the reference point to the candidate point"; + reference_cost = candidate_cost; + accumulated_reference_model_cost_change = + accumulated_candidate_model_cost_change; + } + } } else { ++summary->num_unsuccessful_steps; if (iteration_summary.step_is_valid) {