Fixes for some line search bugs & corner cases.
- Increase precision of numeric values output in error messages to
allow for easier debugging.
- Ensure termination after Wolfe search bracketing phase if bracket
width has been shrunk to below tolerance.
- Cleaned up return value for BracketingPhase(), now false iff
optimisation should stop, true otherwise.
- Fix bug whereby we would mark a step size as satisfying the Wolfe
conditions when it did not due to numerical issues in the cost
function.
- Adding explanation of a subtlety in which a zoom could still be
acceptably invoked with bracket_low.f > bracket_high.f.
- Replacing hard check of a pre-condition of ZoomPhase() with a
conditional return if not satisfied to address issue whereby a
bracket could be incorrectly identified due to inconsistent values
& gradients returned from the cost function.
- Adding missing check for step size validity in line search minimizer.
- Adding ToDebugString() for FunctionSample.
Change-Id: Iad98e635749877f80c079ebad126bf022d82232d
diff --git a/internal/ceres/line_search.cc b/internal/ceres/line_search.cc
index 8323896..2601398 100644
--- a/internal/ceres/line_search.cc
+++ b/internal/ceres/line_search.cc
@@ -29,6 +29,9 @@
// Author: sameeragarwal@google.com (Sameer Agarwal)
#ifndef CERES_NO_LINE_SEARCH_MINIMIZER
+#include <iomanip> // For std::setprecision.
+#include <iostream> // For std::scientific.
+
#include "ceres/line_search.h"
#include "ceres/fpclassify.h"
@@ -41,6 +44,8 @@
namespace ceres {
namespace internal {
namespace {
+// Precision used for floating point values in error message output.
+const int kErrorMessageNumericPrecision = 8;
FunctionSample ValueSample(const double x, const double value) {
FunctionSample sample;
@@ -67,10 +72,7 @@
// Convenience stream operator for pushing FunctionSamples into log messages.
std::ostream& operator<<(std::ostream &os,
const FunctionSample& sample) {
- os << "[x: " << sample.x << ", value: " << sample.value
- << ", gradient: " << sample.gradient << ", value_is_valid: "
- << std::boolalpha << sample.value_is_valid << ", gradient_is_valid: "
- << std::boolalpha << sample.gradient_is_valid << "]";
+ os << sample.ToDebugString();
return os;
}
@@ -170,6 +172,7 @@
// to avoid replicating current.value_is_valid == false
// behaviour in WolfeLineSearch.
CHECK(lowerbound.value_is_valid)
+ << std::scientific << std::setprecision(kErrorMessageNumericPrecision)
<< "Ceres bug: lower-bound sample for interpolation is invalid, "
<< "please contact the developers!, interpolation_type: "
<< LineSearchInterpolationTypeToString(interpolation_type)
@@ -237,20 +240,26 @@
FunctionSample current = ValueAndGradientSample(step_size_estimate, 0.0, 0.0);
current.value_is_valid = false;
- const bool interpolation_uses_gradients =
+ // As the Armijo line search algorithm always uses the initial point, for
+ // which both the function value and derivative are known, when fitting a
+ // minimizing polynomial, we can fit up to a quadratic without requiring the
+ // gradient at the current query point.
+ const bool interpolation_uses_gradient_at_current_sample =
options().interpolation_type == CUBIC;
const double descent_direction_max_norm =
static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
++summary->num_function_evaluations;
- if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
+ if (interpolation_uses_gradient_at_current_sample) {
+ ++summary->num_gradient_evaluations;
+ }
current.value_is_valid =
function->Evaluate(current.x,
¤t.value,
- interpolation_uses_gradients
+ interpolation_uses_gradient_at_current_sample
? ¤t.gradient : NULL);
current.gradient_is_valid =
- interpolation_uses_gradients && current.value_is_valid;
+ interpolation_uses_gradient_at_current_sample && current.value_is_valid;
while (!current.value_is_valid ||
current.value > (initial_cost
+ options().sufficient_decrease
@@ -291,14 +300,16 @@
current.x = step_size;
++summary->num_function_evaluations;
- if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
+ if (interpolation_uses_gradient_at_current_sample) {
+ ++summary->num_gradient_evaluations;
+ }
current.value_is_valid =
function->Evaluate(current.x,
¤t.value,
- interpolation_uses_gradients
+ interpolation_uses_gradient_at_current_sample
? ¤t.gradient : NULL);
current.gradient_is_valid =
- interpolation_uses_gradients && current.value_is_valid;
+ interpolation_uses_gradient_at_current_sample && current.value_is_valid;
}
summary->optimal_step_size = current.x;
@@ -350,28 +361,24 @@
&bracket_low,
&bracket_high,
&do_zoom_search,
- summary) &&
- summary->num_iterations < options().max_num_iterations) {
- // Failed to find either a valid point or a valid bracket, but we did not
- // run out of iterations.
+ summary)) {
+ // Failed to find either a valid point, a valid bracket satisfying the Wolfe
+ // conditions, or even a step size > minimum tolerance satisfying the Armijo
+ // condition.
return;
}
+
if (!do_zoom_search) {
// Either: Bracketing phase already found a point satisfying the strong
// Wolfe conditions, thus no Zoom required.
//
// Or: Bracketing failed to find a valid bracket or a point satisfying the
- // strong Wolfe conditions within max_num_iterations. As this is an
- // 'artificial' constraint, and we would otherwise fail to produce a valid
- // point when ArmijoLineSearch would succeed, we return the lowest point
- // found thus far which satsifies the Armijo condition (but not the Wolfe
- // conditions).
- CHECK(bracket_low.value_is_valid)
- << "Ceres bug: Bracketing produced an invalid bracket_low, please "
- << "contact the developers!, bracket_low: " << bracket_low
- << ", bracket_high: " << bracket_high << ", num_iterations: "
- << summary->num_iterations << ", max_num_iterations: "
- << options().max_num_iterations;
+ // strong Wolfe conditions within max_num_iterations, or whilst searching
+ // shrank the bracket width until it was below our minimum tolerance.
+ // As these are 'artificial' constraints, and we would otherwise fail to
+ // produce a valid point when ArmijoLineSearch would succeed, we return the
+ // point with the lowest cost found thus far which satsifies the Armijo
+ // condition (but not the Wolfe conditions).
summary->optimal_step_size = bracket_low.x;
summary->success = true;
return;
@@ -419,11 +426,22 @@
summary->success = true;
}
-// Returns true iff bracket_low & bracket_high bound a bracket that contains
-// points which satisfy the strong Wolfe conditions. Otherwise, on return false,
-// if we stopped searching due to the 'artificial' condition of reaching
-// max_num_iterations, bracket_low is the step size amongst all those
-// tested, which satisfied the Armijo decrease condition and minimized f().
+// Returns true if either:
+//
+// A termination condition satisfying the (strong) Wolfe bracketing conditions
+// is found:
+//
+// - A valid point, defined as a bracket of zero width [zoom not required].
+// - A valid bracket (of width > tolerance), [zoom required].
+//
+// Or, searching was stopped due to an 'artificial' constraint, i.e. not
+// a condition imposed / required by the underlying algorithm, but instead an
+// engineering / implementation consideration. But a step which exceeds the
+// minimum step size, and satsifies the Armijo condition was still found,
+// and should thus be used [zoom not required].
+//
+// Returns false if no step size > minimum step size was found which
+// satisfies at least the Armijo condition.
bool WolfeLineSearch::BracketingPhase(
const FunctionSample& initial_position,
const double step_size_estimate,
@@ -437,23 +455,28 @@
FunctionSample current = ValueAndGradientSample(step_size_estimate, 0.0, 0.0);
current.value_is_valid = false;
- const bool interpolation_uses_gradients =
- options().interpolation_type == CUBIC;
const double descent_direction_max_norm =
static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
*do_zoom_search = false;
*bracket_low = initial_position;
+ // As we require the gradient to evaluate the Wolfe condition, we always
+ // calculate it together with the value, irrespective of the interpolation
+ // type. As opposed to only calculating the gradient after the Armijo
+ // condition is satisifed, as the computational saving from this approach
+ // would be slight (perhaps even negative due to the extra call). Also,
+ // always calculating the value & gradient together protects against us
+ // reporting invalid solutions if the cost function returns slightly different
+ // function values when evaluated with / without gradients (due to numerical
+ // issues).
++summary->num_function_evaluations;
- if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
+ ++summary->num_gradient_evaluations;
current.value_is_valid =
function->Evaluate(current.x,
¤t.value,
- interpolation_uses_gradients
- ? ¤t.gradient : NULL);
- current.gradient_is_valid =
- interpolation_uses_gradients && current.value_is_valid;
+ ¤t.gradient);
+ current.gradient_is_valid = current.value_is_valid;
while (true) {
++summary->num_iterations;
@@ -473,19 +496,6 @@
break;
}
- // Irrespective of the interpolation type we are using, we now need the
- // gradient at the current point (which satisfies the Armijo condition)
- // in order to check the strong Wolfe conditions.
- if (!interpolation_uses_gradients) {
- ++summary->num_function_evaluations;
- ++summary->num_gradient_evaluations;
- current.value_is_valid =
- function->Evaluate(current.x,
- ¤t.value,
- ¤t.gradient);
- current.gradient_is_valid = current.value_is_valid;
- }
-
if (current.value_is_valid &&
fabs(current.gradient) <=
-options().sufficient_curvature_decrease * initial_position.gradient) {
@@ -507,6 +517,26 @@
*bracket_high = previous;
break;
+ } else if (current.value_is_valid &&
+ fabs(current.x - previous.x) * descent_direction_max_norm
+ < options().min_step_size) {
+ // We have shrunk the search bracket to a width less than our tolerance,
+ // and still not found either a point satisfying the strong Wolfe
+ // conditions, or a valid bracket containing such a point. Stop searching
+ // and set bracket_low to the size size amongst all those tested which
+ // minimizes f() and satisfies the Armijo condition.
+ LOG(WARNING) << "Line search failed: Wolfe bracketing phase shrank "
+ << "bracket width: " << fabs(current.x - previous.x)
+ << ", to < tolerance: " << options().min_step_size
+ << ", with descent_direction_max_norm: "
+ << descent_direction_max_norm << ", and failed to find "
+ << "a point satisfying the strong Wolfe conditions or a "
+ << "bracketing containing such a point. Accepting "
+ << "point found satisfying Armijo condition only, to "
+ << "allow continuation.";
+ *bracket_low = current;
+ break;
+
} else if (summary->num_iterations >= options().max_num_iterations) {
// Check num iterations bound here so that we always evaluate the
// max_num_iterations-th iteration against all conditions, and
@@ -523,7 +553,7 @@
*bracket_low =
current.value_is_valid && current.value < bracket_low->value
? current : *bracket_low;
- return false;
+ break;
}
// Either: f(current) is invalid; or, f(current) is valid, but does not
// satisfy the strong Wolfe conditions itself, or the conditions for
@@ -563,17 +593,22 @@
current.x = step_size;
++summary->num_function_evaluations;
- if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
+ ++summary->num_gradient_evaluations;
current.value_is_valid =
function->Evaluate(current.x,
¤t.value,
- interpolation_uses_gradients
- ? ¤t.gradient : NULL);
- current.gradient_is_valid =
- interpolation_uses_gradients && current.value_is_valid;
+ ¤t.gradient);
+ current.gradient_is_valid = current.value_is_valid;
}
- // Either we have a valid point, defined as a bracket of zero width, in which
- // case no zoom is required, or a valid bracket in which to zoom.
+
+ // Ensure that even if a valid bracket was found, we will only mark a zoom
+ // as required if the bracket's width is greater than our minimum tolerance.
+ if (*do_zoom_search &&
+ fabs(bracket_high->x - bracket_low->x) * descent_direction_max_norm
+ < options().min_step_size) {
+ *do_zoom_search = false;
+ }
+
return true;
}
@@ -589,6 +624,7 @@
Function* function = options().function;
CHECK(bracket_low.value_is_valid && bracket_low.gradient_is_valid)
+ << std::scientific << std::setprecision(kErrorMessageNumericPrecision)
<< "Ceres bug: f_low input to Wolfe Zoom invalid, please contact "
<< "the developers!, initial_position: " << initial_position
<< ", bracket_low: " << bracket_low
@@ -599,22 +635,46 @@
// not have been calculated (if bracket_high.value does not satisfy the
// Armijo sufficient decrease condition and interpolation method does not
// require it).
+ //
+ // We also do not require that: bracket_low.value < bracket_high.value,
+ // although this is typical. This is to deal with the case when
+ // bracket_low = initial_position, bracket_high is the first sample,
+ // and bracket_high does not satisfy the Armijo condition, but still has
+ // bracket_high.value < initial_position.value.
CHECK(bracket_high.value_is_valid)
+ << std::scientific << std::setprecision(kErrorMessageNumericPrecision)
<< "Ceres bug: f_high input to Wolfe Zoom invalid, please "
<< "contact the developers!, initial_position: " << initial_position
<< ", bracket_low: " << bracket_low
<< ", bracket_high: "<< bracket_high;
- CHECK_LT(bracket_low.gradient *
- (bracket_high.x - bracket_low.x), 0.0)
- << "Ceres bug: f_high input to Wolfe Zoom does not satisfy gradient "
- << "condition combined with f_low, please contact the developers!"
- << ", initial_position: " << initial_position
- << ", bracket_low: " << bracket_low
- << ", bracket_high: "<< bracket_high;
+
+ if (bracket_low.gradient * (bracket_high.x - bracket_low.x) >= 0) {
+ // The third condition for a valid initial bracket:
+ //
+ // 3. bracket_high is chosen after bracket_low, s.t.
+ // bracket_low.gradient * (bracket_high.x - bracket_low.x) < 0.
+ //
+ // is not satisfied. As this can happen when the users' cost function
+ // returns inconsistent gradient values relative to the function values,
+ // we do not CHECK_LT(), but we do stop processing and return an invalid
+ // value.
+ summary->error =
+ StringPrintf("Line search failed: Wolfe zoom phase passed a bracket "
+ "which does not satisfy: bracket_low.gradient * "
+ "(bracket_high.x - bracket_low.x) < 0 [%.8e !< 0] "
+ "with initial_position: %s, bracket_low: %s, bracket_high:"
+ " %s, the most likely cause of which is the cost function "
+ "returning inconsistent gradient & function values.",
+ bracket_low.gradient * (bracket_high.x - bracket_low.x),
+ initial_position.ToDebugString().c_str(),
+ bracket_low.ToDebugString().c_str(),
+ bracket_high.ToDebugString().c_str());
+ LOG(WARNING) << summary->error;
+ solution->value_is_valid = false;
+ return false;
+ }
const int num_bracketing_iterations = summary->num_iterations;
- const bool interpolation_uses_gradients =
- options().interpolation_type == CUBIC;
const double descent_direction_max_norm =
static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
@@ -669,15 +729,23 @@
upper_bound_step.x);
// No check on magnitude of step size being too small here as it is
// lower-bounded by the initial bracket start point, which was valid.
+ //
+ // As we require the gradient to evaluate the Wolfe condition, we always
+ // calculate it together with the value, irrespective of the interpolation
+ // type. As opposed to only calculating the gradient after the Armijo
+ // condition is satisifed, as the computational saving from this approach
+ // would be slight (perhaps even negative due to the extra call). Also,
+ // always calculating the value & gradient together protects against us
+ // reporting invalid solutions if the cost function returns slightly
+ // different function values when evaluated with / without gradients (due
+ // to numerical issues).
++summary->num_function_evaluations;
- if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
+ ++summary->num_gradient_evaluations;
solution->value_is_valid =
function->Evaluate(solution->x,
&solution->value,
- interpolation_uses_gradients
- ? &solution->gradient : NULL);
- solution->gradient_is_valid =
- interpolation_uses_gradients && solution->value_is_valid;
+ &solution->gradient);
+ solution->gradient_is_valid = solution->value_is_valid;
if (!solution->value_is_valid) {
summary->error =
StringPrintf("Line search failed: Wolfe Zoom phase found "
@@ -701,28 +769,6 @@
}
// Armijo sufficient decrease satisfied, check strong Wolfe condition.
- if (!interpolation_uses_gradients) {
- // Irrespective of the interpolation type we are using, we now need the
- // gradient at the current point (which satisfies the Armijo condition)
- // in order to check the strong Wolfe conditions.
- ++summary->num_function_evaluations;
- ++summary->num_gradient_evaluations;
- solution->value_is_valid =
- function->Evaluate(solution->x,
- &solution->value,
- &solution->gradient);
- solution->gradient_is_valid = solution->value_is_valid;
- if (!solution->value_is_valid) {
- summary->error =
- StringPrintf("Line search failed: Wolfe Zoom phase found "
- "step_size: %.5e, for which function is invalid, "
- "between low_step: %.5e and high_step: %.5e "
- "at which function is valid.",
- solution->x, bracket_low.x, bracket_high.x);
- LOG(WARNING) << summary->error;
- return false;
- }
- }
if (fabs(solution->gradient) <=
-options().sufficient_curvature_decrease * initial_position.gradient) {
// Found a valid termination point satisfying strong Wolfe conditions.
diff --git a/internal/ceres/line_search_minimizer.cc b/internal/ceres/line_search_minimizer.cc
index b7e96c8..1093bad 100644
--- a/internal/ceres/line_search_minimizer.cc
+++ b/internal/ceres/line_search_minimizer.cc
@@ -314,6 +314,18 @@
current_state.cost,
current_state.directional_derivative,
&line_search_summary);
+ if (!line_search_summary.success) {
+ summary->error =
+ StringPrintf("Numerical failure in line search, failed to find "
+ "a valid step size, (did not run out of iterations) "
+ "using initial_step_size: %.5e, initial_cost: %.5e, "
+ "initial_gradient: %.5e.",
+ initial_step_size, current_state.cost,
+ current_state.directional_derivative);
+ LOG_IF(WARNING, is_not_silent) << summary->error;
+ summary->termination_type = NUMERICAL_FAILURE;
+ break;
+ }
current_state.step_size = line_search_summary.optimal_step_size;
delta = current_state.step_size * current_state.search_direction;
@@ -323,6 +335,13 @@
WallTimeInSeconds() - iteration_start_time;
// TODO(sameeragarwal): Collect stats.
+ //
+ // TODO(sameeragarwal): This call to Plus() directly updates the parameter
+ // vector via the VectorRef x. This is incorrect as we check the
+ // gradient and cost changes to determine if the step is accepted
+ // later, as such we could mutate x with a step that is not
+ // subsequently accepted, thus it is possible that
+ // summary->iterations.end()->x != x at termination.
if (!evaluator->Plus(x.data(), delta.data(), x_plus_delta.data())) {
LOG_IF(WARNING, is_not_silent)
<< "x_plus_delta = Plus(x, delta) failed. ";
diff --git a/internal/ceres/polynomial.cc b/internal/ceres/polynomial.cc
index 3238b89..feec884 100644
--- a/internal/ceres/polynomial.cc
+++ b/internal/ceres/polynomial.cc
@@ -37,6 +37,7 @@
#include "Eigen/Dense"
#include "ceres/internal/port.h"
+#include "ceres/stringprintf.h"
#include "glog/logging.h"
namespace ceres {
@@ -255,6 +256,12 @@
}
}
+string FunctionSample::ToDebugString() const {
+ return StringPrintf("[x: %.8e, value: %.8e, gradient: %.8e, "
+ "value_is_valid: %d, gradient_is_valid: %d]",
+ x, value, gradient, value_is_valid, gradient_is_valid);
+}
+
Vector FindInterpolatingPolynomial(const vector<FunctionSample>& samples) {
const int num_samples = samples.size();
int num_constraints = 0;
diff --git a/internal/ceres/polynomial.h b/internal/ceres/polynomial.h
index 42ffdcb..80ce77e 100644
--- a/internal/ceres/polynomial.h
+++ b/internal/ceres/polynomial.h
@@ -95,6 +95,7 @@
gradient(0.0),
gradient_is_valid(false) {
}
+ string ToDebugString() const;
double x;
double value; // value = f(x)
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index 0445cfb..1420fd8 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -124,7 +124,10 @@
dense_linear_algebra_library_type(EIGEN),
sparse_linear_algebra_library_type(SUITE_SPARSE),
line_search_direction_type(LBFGS),
- line_search_type(ARMIJO) {
+ line_search_type(ARMIJO),
+ line_search_interpolation_type(BISECTION),
+ nonlinear_conjugate_gradient_type(FLETCHER_REEVES),
+ max_lbfgs_rank(-1) {
}
string Solver::Summary::BriefReport() const {
