Add line search to the trust region minimizer.
If enabled, the TrustRegionMinimizer can now further improve the
quality of the step with a Armijo line search.
This is the first step towards adding support for bounds on
variables.
Change-Id: I453b42853cfa6ca4f75812900c13d286a473c2df
diff --git a/internal/ceres/minimizer.h b/internal/ceres/minimizer.h
index ee77726..f9f2b51 100644
--- a/internal/ceres/minimizer.h
+++ b/internal/ceres/minimizer.h
@@ -114,6 +114,7 @@
callbacks = options.callbacks;
inner_iteration_minimizer = NULL;
inner_iteration_tolerance = options.inner_iteration_tolerance;
+ is_constrained = false;
}
int max_num_iterations;
@@ -180,6 +181,9 @@
Minimizer* inner_iteration_minimizer;
double inner_iteration_tolerance;
+
+ // Use a bounds constrained optimization algorithm.
+ bool is_constrained;
};
static bool RunCallbacks(const vector<IterationCallback*> callbacks,
diff --git a/internal/ceres/trust_region_minimizer.cc b/internal/ceres/trust_region_minimizer.cc
index f9f14e0..72e89b9 100644
--- a/internal/ceres/trust_region_minimizer.cc
+++ b/internal/ceres/trust_region_minimizer.cc
@@ -44,6 +44,7 @@
#include "ceres/file.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/scoped_ptr.h"
+#include "ceres/line_search.h"
#include "ceres/linear_least_squares_problems.h"
#include "ceres/sparse_matrix.h"
#include "ceres/stringprintf.h"
@@ -57,6 +58,45 @@
namespace {
// Small constant for various floating point issues.
const double kEpsilon = 1e-12;
+
+LineSearch::Summary DoLineSearch(const Minimizer::Options& options,
+ const Vector& x,
+ const Vector& gradient,
+ const double cost,
+ const Vector& delta,
+ Evaluator* evaluator) {
+ LineSearchFunction line_search_function(evaluator);
+ line_search_function.Init(x, delta);
+
+ LineSearch::Options line_search_options;
+ line_search_options.interpolation_type =
+ options.line_search_interpolation_type;
+ line_search_options.min_step_size = options.min_line_search_step_size;
+ line_search_options.sufficient_decrease =
+ options.line_search_sufficient_function_decrease;
+ line_search_options.max_step_contraction =
+ options.max_line_search_step_contraction;
+ line_search_options.min_step_contraction =
+ options.min_line_search_step_contraction;
+ line_search_options.max_num_iterations =
+ options.max_num_line_search_step_size_iterations;
+ line_search_options.sufficient_curvature_decrease =
+ options.line_search_sufficient_curvature_decrease;
+ line_search_options.max_step_expansion =
+ options.max_line_search_step_expansion;
+ line_search_options.function = &line_search_function;
+
+ string message;
+ scoped_ptr<LineSearch>
+ line_search(CHECK_NOTNULL(
+ LineSearch::Create(ceres::ARMIJO,
+ line_search_options,
+ &message)));
+ LineSearch::Summary summary;
+ line_search->Search(1.0, cost, gradient.dot(delta), &summary);
+ return summary;
+}
+
} // namespace
// Compute a scaling vector that is used to improve the conditioning
@@ -81,24 +121,30 @@
double start_time = WallTimeInSeconds();
double iteration_start_time = start_time;
Init(options);
- const bool is_not_silent = !options.is_silent;
-
- 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 bool is_not_silent = !options.is_silent;
+
+ // If the problem is bounds constrained, then enable the use of a
+ // line search after the trust region step has been computed. This
+ // line search will automatically use a projected the test point
+ // onto the feasible set, there by guaranteeing the feasibility of
+ // the final output.
+ //
+ // TODO(sameeragarwal): Make line search available more generally.
+ const bool use_line_search = options.is_constrained;
+
+ summary->termination_type = NO_CONVERGENCE;
+ summary->num_successful_steps = 0;
+ summary->num_unsuccessful_steps = 0;
+
const int num_parameters = evaluator->NumParameters();
const int num_effective_parameters = evaluator->NumEffectiveParameters();
const int num_residuals = evaluator->NumResiduals();
- VectorRef x_min(parameters, num_parameters);
- Vector x = x_min;
- double x_norm = x.norm();
-
Vector residuals(num_residuals);
Vector trust_region_step(num_effective_parameters);
Vector delta(num_effective_parameters);
@@ -121,6 +167,23 @@
iteration_summary.linear_solver_iterations = 0;
iteration_summary.step_solver_time_in_seconds = 0;
+ VectorRef x_min(parameters, num_parameters);
+ Vector x = x_min;
+ // Project onto the feasible set.
+ if (options.is_constrained) {
+ delta.setZero();
+ if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
+ summary->message = "Unable to project initial point onto the feasible set.";
+ summary->termination_type = FAILURE;
+ LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
+ return;
+ }
+ x_min = x_plus_delta;
+ x = x_plus_delta;
+ }
+
+ double x_norm = x.norm();
+
// Do initial cost and Jacobian evaluation.
double cost = 0.0;
if (!evaluator->Evaluate(x.data(),
@@ -128,9 +191,9 @@
residuals.data(),
gradient.data(),
jacobian)) {
- summary->message = "Terminating: Residual and Jacobian evaluation failed.";
+ summary->message = "Residual and Jacobian evaluation failed.";
summary->termination_type = FAILURE;
- LOG_IF(WARNING, is_not_silent) << summary->message;
+ LOG_IF(WARNING, is_not_silent) << "Terminating: " << summary->message;
return;
}
@@ -305,19 +368,36 @@
// Undo the Jacobian column scaling.
delta = (trust_region_step.array() * scale.array()).matrix();
- double new_cost = numeric_limits<double>::max();
- if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
+ // Try improving the step further by using an ARMIJO line
+ // search.
+ //
+ // TODO(sameeragarwal): What happens to trust region sizing as
+ // it interacts with the line search ?
+ if (use_line_search) {
+ const LineSearch::Summary line_search_summary =
+ DoLineSearch(options, x, gradient, cost, delta, evaluator);
+ if (line_search_summary.success) {
+ delta *= line_search_summary.optimal_step_size;
+ }
+ }
+
+ double new_cost = std::numeric_limits<double>::max();
+ if (evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
+ if (!evaluator->Evaluate(x_plus_delta.data(),
+ &new_cost,
+ NULL,
+ NULL,
+ NULL)) {
+ LOG(WARNING) << "Step failed to evaluate. "
+ << "Treating it as a step with infinite cost";
+ new_cost = numeric_limits<double>::max();
+ }
+ } else {
LOG(WARNING) << "x_plus_delta = Plus(x, delta) failed. "
<< "Treating it as a step with infinite cost";
- } else if (!evaluator->Evaluate(x_plus_delta.data(),
- &new_cost,
- NULL,
- NULL,
- NULL)) {
- LOG(WARNING) << "Step failed to evaluate. "
- << "Treating it as a step with infinite cost";
- new_cost = numeric_limits<double>::max();
- } else {
+ }
+
+ if (new_cost < std::numeric_limits<double>::max()) {
// Check if performing an inner iteration will make it better.
if (inner_iterations_are_enabled) {
++summary->num_inner_iteration_steps;
