Home | History | Annotate | Download | only in ceres
      1 // Ceres Solver - A fast non-linear least squares minimizer
      2 // Copyright 2012 Google Inc. All rights reserved.
      3 // http://code.google.com/p/ceres-solver/
      4 //
      5 // Redistribution and use in source and binary forms, with or without
      6 // modification, are permitted provided that the following conditions are met:
      7 //
      8 // * Redistributions of source code must retain the above copyright notice,
      9 //   this list of conditions and the following disclaimer.
     10 // * Redistributions in binary form must reproduce the above copyright notice,
     11 //   this list of conditions and the following disclaimer in the documentation
     12 //   and/or other materials provided with the distribution.
     13 // * Neither the name of Google Inc. nor the names of its contributors may be
     14 //   used to endorse or promote products derived from this software without
     15 //   specific prior written permission.
     16 //
     17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
     18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
     21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
     22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
     23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
     24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
     25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
     26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
     27 // POSSIBILITY OF SUCH DAMAGE.
     28 //
     29 // Author: sameeragarwal (at) google.com (Sameer Agarwal)
     30 
     31 #ifndef CERES_NO_LINE_SEARCH_MINIMIZER
     32 #include "ceres/line_search.h"
     33 
     34 #include "ceres/fpclassify.h"
     35 #include "ceres/evaluator.h"
     36 #include "ceres/internal/eigen.h"
     37 #include "ceres/polynomial.h"
     38 #include "ceres/stringprintf.h"
     39 #include "glog/logging.h"
     40 
     41 namespace ceres {
     42 namespace internal {
     43 namespace {
     44 
     45 FunctionSample ValueSample(const double x, const double value) {
     46   FunctionSample sample;
     47   sample.x = x;
     48   sample.value = value;
     49   sample.value_is_valid = true;
     50   return sample;
     51 };
     52 
     53 FunctionSample ValueAndGradientSample(const double x,
     54                                       const double value,
     55                                       const double gradient) {
     56   FunctionSample sample;
     57   sample.x = x;
     58   sample.value = value;
     59   sample.gradient = gradient;
     60   sample.value_is_valid = true;
     61   sample.gradient_is_valid = true;
     62   return sample;
     63 };
     64 
     65 }  // namespace
     66 
     67 // Convenience stream operator for pushing FunctionSamples into log messages.
     68 std::ostream& operator<<(std::ostream &os,
     69                          const FunctionSample& sample) {
     70   os << "[x: " << sample.x << ", value: " << sample.value
     71      << ", gradient: " << sample.gradient << ", value_is_valid: "
     72      << std::boolalpha << sample.value_is_valid << ", gradient_is_valid: "
     73      << std::boolalpha << sample.gradient_is_valid << "]";
     74   return os;
     75 }
     76 
     77 LineSearch::LineSearch(const LineSearch::Options& options)
     78     : options_(options) {}
     79 
     80 LineSearch* LineSearch::Create(const LineSearchType line_search_type,
     81                                const LineSearch::Options& options,
     82                                string* error) {
     83   LineSearch* line_search = NULL;
     84   switch (line_search_type) {
     85   case ceres::ARMIJO:
     86     line_search = new ArmijoLineSearch(options);
     87     break;
     88   case ceres::WOLFE:
     89     line_search = new WolfeLineSearch(options);
     90     break;
     91   default:
     92     *error = string("Invalid line search algorithm type: ") +
     93         LineSearchTypeToString(line_search_type) +
     94         string(", unable to create line search.");
     95     return NULL;
     96   }
     97   return line_search;
     98 }
     99 
    100 LineSearchFunction::LineSearchFunction(Evaluator* evaluator)
    101     : evaluator_(evaluator),
    102       position_(evaluator->NumParameters()),
    103       direction_(evaluator->NumEffectiveParameters()),
    104       evaluation_point_(evaluator->NumParameters()),
    105       scaled_direction_(evaluator->NumEffectiveParameters()),
    106       gradient_(evaluator->NumEffectiveParameters()) {
    107 }
    108 
    109 void LineSearchFunction::Init(const Vector& position,
    110                               const Vector& direction) {
    111   position_ = position;
    112   direction_ = direction;
    113 }
    114 
    115 bool LineSearchFunction::Evaluate(double x, double* f, double* g) {
    116   scaled_direction_ = x * direction_;
    117   if (!evaluator_->Plus(position_.data(),
    118                         scaled_direction_.data(),
    119                         evaluation_point_.data())) {
    120     return false;
    121   }
    122 
    123   if (g == NULL) {
    124     return (evaluator_->Evaluate(evaluation_point_.data(),
    125                                   f, NULL, NULL, NULL) &&
    126             IsFinite(*f));
    127   }
    128 
    129   if (!evaluator_->Evaluate(evaluation_point_.data(),
    130                             f,
    131                             NULL,
    132                             gradient_.data(), NULL)) {
    133     return false;
    134   }
    135 
    136   *g = direction_.dot(gradient_);
    137   return IsFinite(*f) && IsFinite(*g);
    138 }
    139 
    140 double LineSearchFunction::DirectionInfinityNorm() const {
    141   return direction_.lpNorm<Eigen::Infinity>();
    142 }
    143 
    144 // Returns step_size \in [min_step_size, max_step_size] which minimizes the
    145 // polynomial of degree defined by interpolation_type which interpolates all
    146 // of the provided samples with valid values.
    147 double LineSearch::InterpolatingPolynomialMinimizingStepSize(
    148     const LineSearchInterpolationType& interpolation_type,
    149     const FunctionSample& lowerbound,
    150     const FunctionSample& previous,
    151     const FunctionSample& current,
    152     const double min_step_size,
    153     const double max_step_size) const {
    154   if (!current.value_is_valid ||
    155       (interpolation_type == BISECTION &&
    156        max_step_size <= current.x)) {
    157     // Either: sample is invalid; or we are using BISECTION and contracting
    158     // the step size.
    159     return min(max(current.x * 0.5, min_step_size), max_step_size);
    160   } else if (interpolation_type == BISECTION) {
    161     CHECK_GT(max_step_size, current.x);
    162     // We are expanding the search (during a Wolfe bracketing phase) using
    163     // BISECTION interpolation.  Using BISECTION when trying to expand is
    164     // strictly speaking an oxymoron, but we define this to mean always taking
    165     // the maximum step size so that the Armijo & Wolfe implementations are
    166     // agnostic to the interpolation type.
    167     return max_step_size;
    168   }
    169   // Only check if lower-bound is valid here, where it is required
    170   // to avoid replicating current.value_is_valid == false
    171   // behaviour in WolfeLineSearch.
    172   CHECK(lowerbound.value_is_valid)
    173       << "Ceres bug: lower-bound sample for interpolation is invalid, "
    174       << "please contact the developers!, interpolation_type: "
    175       << LineSearchInterpolationTypeToString(interpolation_type)
    176       << ", lowerbound: " << lowerbound << ", previous: " << previous
    177       << ", current: " << current;
    178 
    179   // Select step size by interpolating the function and gradient values
    180   // and minimizing the corresponding polynomial.
    181   vector<FunctionSample> samples;
    182   samples.push_back(lowerbound);
    183 
    184   if (interpolation_type == QUADRATIC) {
    185     // Two point interpolation using function values and the
    186     // gradient at the lower bound.
    187     samples.push_back(ValueSample(current.x, current.value));
    188 
    189     if (previous.value_is_valid) {
    190       // Three point interpolation, using function values and the
    191       // gradient at the lower bound.
    192       samples.push_back(ValueSample(previous.x, previous.value));
    193     }
    194   } else if (interpolation_type == CUBIC) {
    195     // Two point interpolation using the function values and the gradients.
    196     samples.push_back(current);
    197 
    198     if (previous.value_is_valid) {
    199       // Three point interpolation using the function values and
    200       // the gradients.
    201       samples.push_back(previous);
    202     }
    203   } else {
    204     LOG(FATAL) << "Ceres bug: No handler for interpolation_type: "
    205                << LineSearchInterpolationTypeToString(interpolation_type)
    206                << ", please contact the developers!";
    207   }
    208 
    209   double step_size = 0.0, unused_min_value = 0.0;
    210   MinimizeInterpolatingPolynomial(samples, min_step_size, max_step_size,
    211                                   &step_size, &unused_min_value);
    212   return step_size;
    213 }
    214 
    215 ArmijoLineSearch::ArmijoLineSearch(const LineSearch::Options& options)
    216     : LineSearch(options) {}
    217 
    218 void ArmijoLineSearch::Search(const double step_size_estimate,
    219                               const double initial_cost,
    220                               const double initial_gradient,
    221                               Summary* summary) {
    222   *CHECK_NOTNULL(summary) = LineSearch::Summary();
    223   CHECK_GE(step_size_estimate, 0.0);
    224   CHECK_GT(options().sufficient_decrease, 0.0);
    225   CHECK_LT(options().sufficient_decrease, 1.0);
    226   CHECK_GT(options().max_num_iterations, 0);
    227   Function* function = options().function;
    228 
    229   // Note initial_cost & initial_gradient are evaluated at step_size = 0,
    230   // not step_size_estimate, which is our starting guess.
    231   const FunctionSample initial_position =
    232       ValueAndGradientSample(0.0, initial_cost, initial_gradient);
    233 
    234   FunctionSample previous = ValueAndGradientSample(0.0, 0.0, 0.0);
    235   previous.value_is_valid = false;
    236 
    237   FunctionSample current = ValueAndGradientSample(step_size_estimate, 0.0, 0.0);
    238   current.value_is_valid = false;
    239 
    240   const bool interpolation_uses_gradients =
    241       options().interpolation_type == CUBIC;
    242   const double descent_direction_max_norm =
    243       static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
    244 
    245   ++summary->num_function_evaluations;
    246   if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
    247   current.value_is_valid =
    248       function->Evaluate(current.x,
    249                          &current.value,
    250                          interpolation_uses_gradients
    251                          ? &current.gradient : NULL);
    252   current.gradient_is_valid =
    253       interpolation_uses_gradients && current.value_is_valid;
    254   while (!current.value_is_valid ||
    255          current.value > (initial_cost
    256                           + options().sufficient_decrease
    257                           * initial_gradient
    258                           * current.x)) {
    259     // If current.value_is_valid is false, we treat it as if the cost at that
    260     // point is not large enough to satisfy the sufficient decrease condition.
    261     ++summary->num_iterations;
    262     if (summary->num_iterations >= options().max_num_iterations) {
    263       summary->error =
    264           StringPrintf("Line search failed: Armijo failed to find a point "
    265                        "satisfying the sufficient decrease condition within "
    266                        "specified max_num_iterations: %d.",
    267                        options().max_num_iterations);
    268       LOG(WARNING) << summary->error;
    269       return;
    270     }
    271 
    272     const double step_size =
    273         this->InterpolatingPolynomialMinimizingStepSize(
    274             options().interpolation_type,
    275             initial_position,
    276             previous,
    277             current,
    278             (options().max_step_contraction * current.x),
    279             (options().min_step_contraction * current.x));
    280 
    281     if (step_size * descent_direction_max_norm < options().min_step_size) {
    282       summary->error =
    283           StringPrintf("Line search failed: step_size too small: %.5e "
    284                        "with descent_direction_max_norm: %.5e.", step_size,
    285                        descent_direction_max_norm);
    286       LOG(WARNING) << summary->error;
    287       return;
    288     }
    289 
    290     previous = current;
    291     current.x = step_size;
    292 
    293     ++summary->num_function_evaluations;
    294     if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
    295     current.value_is_valid =
    296       function->Evaluate(current.x,
    297                          &current.value,
    298                          interpolation_uses_gradients
    299                          ? &current.gradient : NULL);
    300     current.gradient_is_valid =
    301         interpolation_uses_gradients && current.value_is_valid;
    302   }
    303 
    304   summary->optimal_step_size = current.x;
    305   summary->success = true;
    306 }
    307 
    308 WolfeLineSearch::WolfeLineSearch(const LineSearch::Options& options)
    309     : LineSearch(options) {}
    310 
    311 void WolfeLineSearch::Search(const double step_size_estimate,
    312                              const double initial_cost,
    313                              const double initial_gradient,
    314                              Summary* summary) {
    315   *CHECK_NOTNULL(summary) = LineSearch::Summary();
    316   // All parameters should have been validated by the Solver, but as
    317   // invalid values would produce crazy nonsense, hard check them here.
    318   CHECK_GE(step_size_estimate, 0.0);
    319   CHECK_GT(options().sufficient_decrease, 0.0);
    320   CHECK_GT(options().sufficient_curvature_decrease,
    321            options().sufficient_decrease);
    322   CHECK_LT(options().sufficient_curvature_decrease, 1.0);
    323   CHECK_GT(options().max_step_expansion, 1.0);
    324 
    325   // Note initial_cost & initial_gradient are evaluated at step_size = 0,
    326   // not step_size_estimate, which is our starting guess.
    327   const FunctionSample initial_position =
    328       ValueAndGradientSample(0.0, initial_cost, initial_gradient);
    329 
    330   bool do_zoom_search = false;
    331   // Important: The high/low in bracket_high & bracket_low refer to their
    332   // _function_ values, not their step sizes i.e. it is _not_ required that
    333   // bracket_low.x < bracket_high.x.
    334   FunctionSample solution, bracket_low, bracket_high;
    335 
    336   // Wolfe bracketing phase: Increases step_size until either it finds a point
    337   // that satisfies the (strong) Wolfe conditions, or an interval that brackets
    338   // step sizes which satisfy the conditions.  From Nocedal & Wright [1] p61 the
    339   // interval: (step_size_{k-1}, step_size_{k}) contains step lengths satisfying
    340   // the strong Wolfe conditions if one of the following conditions are met:
    341   //
    342   //   1. step_size_{k} violates the sufficient decrease (Armijo) condition.
    343   //   2. f(step_size_{k}) >= f(step_size_{k-1}).
    344   //   3. f'(step_size_{k}) >= 0.
    345   //
    346   // Caveat: If f(step_size_{k}) is invalid, then step_size is reduced, ignoring
    347   // this special case, step_size monotonically increases during bracketing.
    348   if (!this->BracketingPhase(initial_position,
    349                              step_size_estimate,
    350                              &bracket_low,
    351                              &bracket_high,
    352                              &do_zoom_search,
    353                              summary) &&
    354       summary->num_iterations < options().max_num_iterations) {
    355     // Failed to find either a valid point or a valid bracket, but we did not
    356     // run out of iterations.
    357     return;
    358   }
    359   if (!do_zoom_search) {
    360     // Either: Bracketing phase already found a point satisfying the strong
    361     // Wolfe conditions, thus no Zoom required.
    362     //
    363     // Or: Bracketing failed to find a valid bracket or a point satisfying the
    364     // strong Wolfe conditions within max_num_iterations.  As this is an
    365     // 'artificial' constraint, and we would otherwise fail to produce a valid
    366     // point when ArmijoLineSearch would succeed, we return the lowest point
    367     // found thus far which satsifies the Armijo condition (but not the Wolfe
    368     // conditions).
    369     CHECK(bracket_low.value_is_valid)
    370         << "Ceres bug: Bracketing produced an invalid bracket_low, please "
    371         << "contact the developers!, bracket_low: " << bracket_low
    372         << ", bracket_high: " << bracket_high << ", num_iterations: "
    373         << summary->num_iterations << ", max_num_iterations: "
    374         << options().max_num_iterations;
    375     summary->optimal_step_size = bracket_low.x;
    376     summary->success = true;
    377     return;
    378   }
    379 
    380   // Wolfe Zoom phase: Called when the Bracketing phase finds an interval of
    381   // non-zero, finite width that should bracket step sizes which satisfy the
    382   // (strong) Wolfe conditions (before finding a step size that satisfies the
    383   // conditions).  Zoom successively decreases the size of the interval until a
    384   // step size which satisfies the Wolfe conditions is found.  The interval is
    385   // defined by bracket_low & bracket_high, which satisfy:
    386   //
    387   //   1. The interval bounded by step sizes: bracket_low.x & bracket_high.x
    388   //      contains step sizes that satsify the strong Wolfe conditions.
    389   //   2. bracket_low.x is of all the step sizes evaluated *which satisifed the
    390   //      Armijo sufficient decrease condition*, the one which generated the
    391   //      smallest function value, i.e. bracket_low.value <
    392   //      f(all other steps satisfying Armijo).
    393   //        - Note that this does _not_ (necessarily) mean that initially
    394   //          bracket_low.value < bracket_high.value (although this is typical)
    395   //          e.g. when bracket_low = initial_position, and bracket_high is the
    396   //          first sample, and which does not satisfy the Armijo condition,
    397   //          but still has bracket_high.value < initial_position.value.
    398   //   3. bracket_high is chosen after bracket_low, s.t.
    399   //      bracket_low.gradient * (bracket_high.x - bracket_low.x) < 0.
    400   if (!this->ZoomPhase(initial_position,
    401                        bracket_low,
    402                        bracket_high,
    403                        &solution,
    404                        summary) && !solution.value_is_valid) {
    405     // Failed to find a valid point (given the specified decrease parameters)
    406     // within the specified bracket.
    407     return;
    408   }
    409   // Ensure that if we ran out of iterations whilst zooming the bracket, or
    410   // shrank the bracket width to < tolerance and failed to find a point which
    411   // satisfies the strong Wolfe curvature condition, that we return the point
    412   // amongst those found thus far, which minimizes f() and satisfies the Armijo
    413   // condition.
    414   solution =
    415       solution.value_is_valid && solution.value <= bracket_low.value
    416       ? solution : bracket_low;
    417 
    418   summary->optimal_step_size = solution.x;
    419   summary->success = true;
    420 }
    421 
    422 // Returns true iff bracket_low & bracket_high bound a bracket that contains
    423 // points which satisfy the strong Wolfe conditions. Otherwise, on return false,
    424 // if we stopped searching due to the 'artificial' condition of reaching
    425 // max_num_iterations, bracket_low is the step size amongst all those
    426 // tested, which satisfied the Armijo decrease condition and minimized f().
    427 bool WolfeLineSearch::BracketingPhase(
    428     const FunctionSample& initial_position,
    429     const double step_size_estimate,
    430     FunctionSample* bracket_low,
    431     FunctionSample* bracket_high,
    432     bool* do_zoom_search,
    433     Summary* summary) {
    434   Function* function = options().function;
    435 
    436   FunctionSample previous = initial_position;
    437   FunctionSample current = ValueAndGradientSample(step_size_estimate, 0.0, 0.0);
    438   current.value_is_valid = false;
    439 
    440   const bool interpolation_uses_gradients =
    441       options().interpolation_type == CUBIC;
    442   const double descent_direction_max_norm =
    443       static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
    444 
    445   *do_zoom_search = false;
    446   *bracket_low = initial_position;
    447 
    448   ++summary->num_function_evaluations;
    449   if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
    450   current.value_is_valid =
    451       function->Evaluate(current.x,
    452                          &current.value,
    453                          interpolation_uses_gradients
    454                          ? &current.gradient : NULL);
    455   current.gradient_is_valid =
    456       interpolation_uses_gradients && current.value_is_valid;
    457 
    458   while (true) {
    459     ++summary->num_iterations;
    460 
    461     if (current.value_is_valid &&
    462         (current.value > (initial_position.value
    463                           + options().sufficient_decrease
    464                           * initial_position.gradient
    465                           * current.x) ||
    466          (previous.value_is_valid && current.value > previous.value))) {
    467       // Bracket found: current step size violates Armijo sufficient decrease
    468       // condition, or has stepped past an inflection point of f() relative to
    469       // previous step size.
    470       *do_zoom_search = true;
    471       *bracket_low = previous;
    472       *bracket_high = current;
    473       break;
    474     }
    475 
    476     // Irrespective of the interpolation type we are using, we now need the
    477     // gradient at the current point (which satisfies the Armijo condition)
    478     // in order to check the strong Wolfe conditions.
    479     if (!interpolation_uses_gradients) {
    480       ++summary->num_function_evaluations;
    481       ++summary->num_gradient_evaluations;
    482       current.value_is_valid =
    483           function->Evaluate(current.x,
    484                              &current.value,
    485                              &current.gradient);
    486       current.gradient_is_valid = current.value_is_valid;
    487     }
    488 
    489     if (current.value_is_valid &&
    490         fabs(current.gradient) <=
    491         -options().sufficient_curvature_decrease * initial_position.gradient) {
    492       // Current step size satisfies the strong Wolfe conditions, and is thus a
    493       // valid termination point, therefore a Zoom not required.
    494       *bracket_low = current;
    495       *bracket_high = current;
    496       break;
    497 
    498     } else if (current.value_is_valid && current.gradient >= 0) {
    499       // Bracket found: current step size has stepped past an inflection point
    500       // of f(), but Armijo sufficient decrease is still satisfied and
    501       // f(current) is our best minimum thus far.  Remember step size
    502       // monotonically increases, thus previous_step_size < current_step_size
    503       // even though f(previous) > f(current).
    504       *do_zoom_search = true;
    505       // Note inverse ordering from first bracket case.
    506       *bracket_low = current;
    507       *bracket_high = previous;
    508       break;
    509 
    510     } else if (summary->num_iterations >= options().max_num_iterations) {
    511       // Check num iterations bound here so that we always evaluate the
    512       // max_num_iterations-th iteration against all conditions, and
    513       // then perform no additional (unused) evaluations.
    514       summary->error =
    515           StringPrintf("Line search failed: Wolfe bracketing phase failed to "
    516                        "find a point satisfying strong Wolfe conditions, or a "
    517                        "bracket containing such a point within specified "
    518                        "max_num_iterations: %d", options().max_num_iterations);
    519       LOG(WARNING) << summary->error;
    520       // Ensure that bracket_low is always set to the step size amongst all
    521       // those tested which minimizes f() and satisfies the Armijo condition
    522       // when we terminate due to the 'artificial' max_num_iterations condition.
    523       *bracket_low =
    524           current.value_is_valid && current.value < bracket_low->value
    525           ? current : *bracket_low;
    526       return false;
    527     }
    528     // Either: f(current) is invalid; or, f(current) is valid, but does not
    529     // satisfy the strong Wolfe conditions itself, or the conditions for
    530     // being a boundary of a bracket.
    531 
    532     // If f(current) is valid, (but meets no criteria) expand the search by
    533     // increasing the step size.
    534     const double max_step_size =
    535         current.value_is_valid
    536         ? (current.x * options().max_step_expansion) : current.x;
    537 
    538     // We are performing 2-point interpolation only here, but the API of
    539     // InterpolatingPolynomialMinimizingStepSize() allows for up to
    540     // 3-point interpolation, so pad call with a sample with an invalid
    541     // value that will therefore be ignored.
    542     const FunctionSample unused_previous;
    543     DCHECK(!unused_previous.value_is_valid);
    544     // Contracts step size if f(current) is not valid.
    545     const double step_size =
    546         this->InterpolatingPolynomialMinimizingStepSize(
    547             options().interpolation_type,
    548             previous,
    549             unused_previous,
    550             current,
    551             previous.x,
    552             max_step_size);
    553     if (step_size * descent_direction_max_norm < options().min_step_size) {
    554       summary->error =
    555           StringPrintf("Line search failed: step_size too small: %.5e "
    556                        "with descent_direction_max_norm: %.5e", step_size,
    557                        descent_direction_max_norm);
    558       LOG(WARNING) << summary->error;
    559       return false;
    560     }
    561 
    562     previous = current.value_is_valid ? current : previous;
    563     current.x = step_size;
    564 
    565     ++summary->num_function_evaluations;
    566     if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
    567     current.value_is_valid =
    568         function->Evaluate(current.x,
    569                            &current.value,
    570                            interpolation_uses_gradients
    571                            ? &current.gradient : NULL);
    572     current.gradient_is_valid =
    573         interpolation_uses_gradients && current.value_is_valid;
    574   }
    575   // Either we have a valid point, defined as a bracket of zero width, in which
    576   // case no zoom is required, or a valid bracket in which to zoom.
    577   return true;
    578 }
    579 
    580 // Returns true iff solution satisfies the strong Wolfe conditions. Otherwise,
    581 // on return false, if we stopped searching due to the 'artificial' condition of
    582 // reaching max_num_iterations, solution is the step size amongst all those
    583 // tested, which satisfied the Armijo decrease condition and minimized f().
    584 bool WolfeLineSearch::ZoomPhase(const FunctionSample& initial_position,
    585                                 FunctionSample bracket_low,
    586                                 FunctionSample bracket_high,
    587                                 FunctionSample* solution,
    588                                 Summary* summary) {
    589   Function* function = options().function;
    590 
    591   CHECK(bracket_low.value_is_valid && bracket_low.gradient_is_valid)
    592       << "Ceres bug: f_low input to Wolfe Zoom invalid, please contact "
    593       << "the developers!, initial_position: " << initial_position
    594       << ", bracket_low: " << bracket_low
    595       << ", bracket_high: "<< bracket_high;
    596   // We do not require bracket_high.gradient_is_valid as the gradient condition
    597   // for a valid bracket is only dependent upon bracket_low.gradient, and
    598   // in order to minimize jacobian evaluations, bracket_high.gradient may
    599   // not have been calculated (if bracket_high.value does not satisfy the
    600   // Armijo sufficient decrease condition and interpolation method does not
    601   // require it).
    602   CHECK(bracket_high.value_is_valid)
    603       << "Ceres bug: f_high input to Wolfe Zoom invalid, please "
    604       << "contact the developers!, initial_position: " << initial_position
    605       << ", bracket_low: " << bracket_low
    606       << ", bracket_high: "<< bracket_high;
    607   CHECK_LT(bracket_low.gradient *
    608            (bracket_high.x - bracket_low.x), 0.0)
    609       << "Ceres bug: f_high input to Wolfe Zoom does not satisfy gradient "
    610       << "condition combined with f_low, please contact the developers!"
    611       << ", initial_position: " << initial_position
    612       << ", bracket_low: " << bracket_low
    613       << ", bracket_high: "<< bracket_high;
    614 
    615   const int num_bracketing_iterations = summary->num_iterations;
    616   const bool interpolation_uses_gradients =
    617       options().interpolation_type == CUBIC;
    618   const double descent_direction_max_norm =
    619       static_cast<const LineSearchFunction*>(function)->DirectionInfinityNorm();
    620 
    621   while (true) {
    622     // Set solution to bracket_low, as it is our best step size (smallest f())
    623     // found thus far and satisfies the Armijo condition, even though it does
    624     // not satisfy the Wolfe condition.
    625     *solution = bracket_low;
    626     if (summary->num_iterations >= options().max_num_iterations) {
    627       summary->error =
    628           StringPrintf("Line search failed: Wolfe zoom phase failed to "
    629                        "find a point satisfying strong Wolfe conditions "
    630                        "within specified max_num_iterations: %d, "
    631                        "(num iterations taken for bracketing: %d).",
    632                        options().max_num_iterations, num_bracketing_iterations);
    633       LOG(WARNING) << summary->error;
    634       return false;
    635     }
    636     if (fabs(bracket_high.x - bracket_low.x) * descent_direction_max_norm
    637         < options().min_step_size) {
    638       // Bracket width has been reduced below tolerance, and no point satisfying
    639       // the strong Wolfe conditions has been found.
    640       summary->error =
    641           StringPrintf("Line search failed: Wolfe zoom bracket width: %.5e "
    642                        "too small with descent_direction_max_norm: %.5e.",
    643                        fabs(bracket_high.x - bracket_low.x),
    644                        descent_direction_max_norm);
    645       LOG(WARNING) << summary->error;
    646       return false;
    647     }
    648 
    649     ++summary->num_iterations;
    650     // Polynomial interpolation requires inputs ordered according to step size,
    651     // not f(step size).
    652     const FunctionSample& lower_bound_step =
    653         bracket_low.x < bracket_high.x ? bracket_low : bracket_high;
    654     const FunctionSample& upper_bound_step =
    655         bracket_low.x < bracket_high.x ? bracket_high : bracket_low;
    656     // We are performing 2-point interpolation only here, but the API of
    657     // InterpolatingPolynomialMinimizingStepSize() allows for up to
    658     // 3-point interpolation, so pad call with a sample with an invalid
    659     // value that will therefore be ignored.
    660     const FunctionSample unused_previous;
    661     DCHECK(!unused_previous.value_is_valid);
    662     solution->x =
    663         this->InterpolatingPolynomialMinimizingStepSize(
    664             options().interpolation_type,
    665             lower_bound_step,
    666             unused_previous,
    667             upper_bound_step,
    668             lower_bound_step.x,
    669             upper_bound_step.x);
    670     // No check on magnitude of step size being too small here as it is
    671     // lower-bounded by the initial bracket start point, which was valid.
    672     ++summary->num_function_evaluations;
    673     if (interpolation_uses_gradients) { ++summary->num_gradient_evaluations; }
    674     solution->value_is_valid =
    675         function->Evaluate(solution->x,
    676                            &solution->value,
    677                            interpolation_uses_gradients
    678                            ? &solution->gradient : NULL);
    679     solution->gradient_is_valid =
    680         interpolation_uses_gradients && solution->value_is_valid;
    681     if (!solution->value_is_valid) {
    682       summary->error =
    683           StringPrintf("Line search failed: Wolfe Zoom phase found "
    684                        "step_size: %.5e, for which function is invalid, "
    685                        "between low_step: %.5e and high_step: %.5e "
    686                        "at which function is valid.",
    687                        solution->x, bracket_low.x, bracket_high.x);
    688       LOG(WARNING) << summary->error;
    689       return false;
    690     }
    691 
    692     if ((solution->value > (initial_position.value
    693                             + options().sufficient_decrease
    694                             * initial_position.gradient
    695                             * solution->x)) ||
    696         (solution->value >= bracket_low.value)) {
    697       // Armijo sufficient decrease not satisfied, or not better
    698       // than current lowest sample, use as new upper bound.
    699       bracket_high = *solution;
    700       continue;
    701     }
    702 
    703     // Armijo sufficient decrease satisfied, check strong Wolfe condition.
    704     if (!interpolation_uses_gradients) {
    705       // Irrespective of the interpolation type we are using, we now need the
    706       // gradient at the current point (which satisfies the Armijo condition)
    707       // in order to check the strong Wolfe conditions.
    708       ++summary->num_function_evaluations;
    709       ++summary->num_gradient_evaluations;
    710       solution->value_is_valid =
    711           function->Evaluate(solution->x,
    712                              &solution->value,
    713                              &solution->gradient);
    714       solution->gradient_is_valid = solution->value_is_valid;
    715       if (!solution->value_is_valid) {
    716         summary->error =
    717             StringPrintf("Line search failed: Wolfe Zoom phase found "
    718                          "step_size: %.5e, for which function is invalid, "
    719                          "between low_step: %.5e and high_step: %.5e "
    720                          "at which function is valid.",
    721                          solution->x, bracket_low.x, bracket_high.x);
    722         LOG(WARNING) << summary->error;
    723         return false;
    724       }
    725     }
    726     if (fabs(solution->gradient) <=
    727         -options().sufficient_curvature_decrease * initial_position.gradient) {
    728       // Found a valid termination point satisfying strong Wolfe conditions.
    729       break;
    730 
    731     } else if (solution->gradient * (bracket_high.x - bracket_low.x) >= 0) {
    732       bracket_high = bracket_low;
    733     }
    734 
    735     bracket_low = *solution;
    736   }
    737   // Solution contains a valid point which satisfies the strong Wolfe
    738   // conditions.
    739   return true;
    740 }
    741 
    742 }  // namespace internal
    743 }  // namespace ceres
    744 
    745 #endif  // CERES_NO_LINE_SEARCH_MINIMIZER
    746