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@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
41namespace ceres {
42namespace internal {
43namespace {
44
45FunctionSample 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
53FunctionSample 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.
68std::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
77LineSearch::LineSearch(const LineSearch::Options& options)
78    : options_(options) {}
79
80LineSearch* 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
100LineSearchFunction::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
109void LineSearchFunction::Init(const Vector& position,
110                              const Vector& direction) {
111  position_ = position;
112  direction_ = direction;
113}
114
115bool 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
140double 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.
147double 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
215ArmijoLineSearch::ArmijoLineSearch(const LineSearch::Options& options)
216    : LineSearch(options) {}
217
218void 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
308WolfeLineSearch::WolfeLineSearch(const LineSearch::Options& options)
309    : LineSearch(options) {}
310
311void 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().
427bool 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().
584bool 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