Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
line_minimization.h
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 //---------------------------------------------------------------
15 
16 #ifndef dealii_line_minimization_h
17 #define dealii_line_minimization_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/logstream.h>
23 #include <deal.II/base/numbers.h>
24 #include <deal.II/base/utilities.h>
25 
26 #include <deal.II/numerics/history.h>
27 
28 #include <boost/optional.hpp>
29 
30 #include <errno.h>
31 #include <sys/stat.h>
32 
33 #include <fstream>
34 #include <string>
35 
36 
37 DEAL_II_NAMESPACE_OPEN
38 
43 {
52  template <typename NumberType>
53  boost::optional<NumberType>
54  quadratic_fit(const NumberType x_low,
55  const NumberType f_low,
56  const NumberType g_low,
57  const NumberType x_hi,
58  const NumberType f_hi);
59 
69  template <typename NumberType>
70  boost::optional<NumberType>
71  cubic_fit(const NumberType x_low,
72  const NumberType f_low,
73  const NumberType g_low,
74  const NumberType x_hi,
75  const NumberType f_hi,
76  const NumberType g_hi);
77 
85  template <typename NumberType>
86  boost::optional<NumberType>
87  cubic_fit_three_points(const NumberType x_low,
88  const NumberType f_low,
89  const NumberType g_low,
90  const NumberType x_hi,
91  const NumberType f_hi,
92  const NumberType x_rec,
93  const NumberType f_rec);
94 
106  template <typename NumberType>
107  NumberType
108  poly_fit(const NumberType x_low,
109  const NumberType f_low,
110  const NumberType g_low,
111  const NumberType x_hi,
112  const NumberType f_hi,
113  const NumberType g_hi,
114  const FiniteSizeHistory<NumberType> & x_rec,
115  const FiniteSizeHistory<NumberType> & f_rec,
116  const FiniteSizeHistory<NumberType> & g_rec,
117  const std::pair<NumberType, NumberType> bounds);
118 
123  template <typename NumberType>
124  NumberType
125  poly_fit_three_points(const NumberType x_low,
126  const NumberType f_low,
127  const NumberType g_low,
128  const NumberType x_hi,
129  const NumberType f_hi,
130  const NumberType g_hi,
131  const FiniteSizeHistory<NumberType> & x_rec,
132  const FiniteSizeHistory<NumberType> & f_rec,
133  const FiniteSizeHistory<NumberType> & g_rec,
134  const std::pair<NumberType, NumberType> bounds);
135 
136 
323  template <typename NumberType>
324  std::pair<NumberType, unsigned int>
325  line_search(
326  const std::function<std::pair<NumberType, NumberType>(const NumberType x)>
327  & func,
328  const NumberType f0,
329  const NumberType g0,
330  const std::function<
331  NumberType(const NumberType x_low,
332  const NumberType f_low,
333  const NumberType g_low,
334  const NumberType x_hi,
335  const NumberType f_hi,
336  const NumberType g_hi,
337  const FiniteSizeHistory<NumberType> & x_rec,
338  const FiniteSizeHistory<NumberType> & f_rec,
339  const FiniteSizeHistory<NumberType> & g_rec,
340  const std::pair<NumberType, NumberType> bounds)> &interpolate,
341  const NumberType a1,
342  const NumberType eta = 0.9,
343  const NumberType mu = 0.01,
344  const NumberType a_max = std::numeric_limits<NumberType>::max(),
345  const unsigned int max_evaluations = 20,
346  const bool debug_output = false);
347 
348 
349  // ------------------- inline and template functions ----------------
350 
351 
352 #ifndef DOXYGEN
353 
354 
355  template <typename NumberType>
356  boost::optional<NumberType>
357  quadratic_fit(const NumberType x1,
358  const NumberType f1,
359  const NumberType g1,
360  const NumberType x2,
361  const NumberType f2)
362  {
363  Assert(x1 != x2, ExcMessage("Point are the same"));
364  const NumberType denom = (2. * g1 * x2 - 2. * g1 * x1 - 2. * f2 + 2. * f1);
365  if (denom == 0)
366  return boost::none;
367  else
368  return (g1 * (x2 * x2 - x1 * x1) + 2. * (f1 - f2) * x1) / denom;
369  }
370 
371 
372 
373  template <typename NumberType>
374  boost::optional<NumberType>
375  cubic_fit(const NumberType x1,
376  const NumberType f1,
377  const NumberType g1,
378  const NumberType x2,
379  const NumberType f2,
380  const NumberType g2)
381  {
382  Assert(x1 != x2, ExcMessage("Points are the same"));
383  const NumberType beta1 = g1 + g2 - 3. * (f1 - f2) / (x1 - x2);
384  const NumberType s = beta1 * beta1 - g1 * g2;
385  if (s < 0)
386  return boost::none;
387 
388  const NumberType beta2 = std::sqrt(s);
389  const NumberType denom =
390  x1 < x2 ? g2 - g1 + 2. * beta2 : g1 - g2 + 2. * beta2;
391  if (denom == 0.)
392  return boost::none;
393 
394  return x1 < x2 ? x2 - (x2 - x1) * (g2 + beta2 - beta1) / denom :
395  x1 - (x1 - x2) * (g1 + beta2 - beta1) / denom;
396  }
397 
398 
399 
400  template <typename NumberType>
401  boost::optional<NumberType>
402  cubic_fit_three_points(const NumberType x1,
403  const NumberType f1,
404  const NumberType g1,
405  const NumberType x2,
406  const NumberType f2,
407  const NumberType x3,
408  const NumberType f3)
409  {
410  Assert(x1 != x2, ExcMessage("Points are the same"));
411  Assert(x1 != x3, ExcMessage("Points are the same"));
412  // f(x) = A *(x-x1)^3 + B*(x-x1)^2 + C*(x-x1) + D
413  // =>
414  // D = f1
415  // C = g1
416 
417  // the rest is a system of 2 equations:
418 
419  const NumberType x2_shift = x2 - x1;
420  const NumberType x3_shift = x3 - x1;
421  const NumberType r1 = f2 - f1 - g1 * x2_shift;
422  const NumberType r2 = f3 - f1 - g1 * x3_shift;
423  const NumberType denom =
424  std::pow(x2_shift * x3_shift, 2) * (x2_shift - x3_shift);
425  if (denom == 0.)
426  return boost::none;
427 
428  const NumberType A =
429  (r1 * std::pow(x3_shift, 2) - r2 * std::pow(x2_shift, 2)) / denom;
430  const NumberType B =
431  (r2 * std::pow(x2_shift, 3) - r1 * std::pow(x3_shift, 3)) / denom;
432  const NumberType &C = g1;
433 
434  // now get the minimizer:
435  const NumberType radical = B * B - A * C * 3;
436  if (radical < 0)
437  return boost::none;
438 
439  return x1 + (-B + std::sqrt(radical)) / (A * 3);
440  }
441 
442 
443 
444  template <typename NumberType>
445  NumberType
446  poly_fit(const NumberType x1,
447  const NumberType f1,
448  const NumberType g1,
449  const NumberType x2,
450  const NumberType f2,
451  const NumberType g2,
455  const std::pair<NumberType, NumberType> bounds)
456  {
457  Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds"));
458 
459  // Similar to scipy implementation but we fit based on two points
460  // with their gradients and do bisection on bounds.
461  // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563
462 
463  // First try cubic interpolation
464  boost::optional<NumberType> res = cubic_fit(x1, f1, g1, x2, f2, g2);
465  if (res && *res >= bounds.first && *res <= bounds.second)
466  return *res;
467 
468  // cubic either fails or outside of safe region, do quadratic:
469  res = quadratic_fit(x1, f1, g1, x2, f2);
470  if (res && *res >= bounds.first && *res <= bounds.second)
471  return *res;
472 
473  // quadratic either failed or outside of safe region. Do bisection
474  // on safe region
475  return (bounds.first + bounds.second) * 0.5;
476  }
477 
478 
479 
480  template <typename NumberType>
481  NumberType
482  poly_fit_three_points(const NumberType x1,
483  const NumberType f1,
484  const NumberType g1,
485  const NumberType x2,
486  const NumberType f2,
487  const NumberType g2,
488  const FiniteSizeHistory<NumberType> &x_rec,
489  const FiniteSizeHistory<NumberType> &f_rec,
490  const FiniteSizeHistory<NumberType> & /*g_rec*/,
491  const std::pair<NumberType, NumberType> bounds)
492  {
493  Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds"));
494  AssertDimension(x_rec.size(), f_rec.size());
495 
496  // Same as scipy implementation where cubic fit is using 3 points
497  // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563
498 
499  // First try cubic interpolation after first iteration
500  boost::optional<NumberType> res =
501  x_rec.size() > 0 ?
502  cubic_fit_three_points(x1, f1, g1, x2, f2, x_rec[0], f_rec[0]) :
503  boost::none;
504  if (res && *res >= bounds.first && *res <= bounds.second)
505  return *res;
506 
507  // cubic either fails or outside of safe region, do quadratic:
508  res = quadratic_fit(x1, f1, g1, x2, f2);
509  if (res && *res >= bounds.first && *res <= bounds.second)
510  return *res;
511 
512  // quadratic either failed or outside of safe region. Do bisection
513  // on safe region
514  return (bounds.first + bounds.second) * 0.5;
515  }
516 
517 
518 
519  template <typename NumberType>
520  std::pair<NumberType, unsigned int>
521  line_search(
522  const std::function<std::pair<NumberType, NumberType>(const NumberType x)>
523  & func,
524  const NumberType f0,
525  const NumberType g0,
526  const std::function<
527  NumberType(const NumberType x_low,
528  const NumberType f_low,
529  const NumberType g_low,
530  const NumberType x_hi,
531  const NumberType f_hi,
532  const NumberType g_hi,
533  const FiniteSizeHistory<NumberType> & x_rec,
534  const FiniteSizeHistory<NumberType> & f_rec,
535  const FiniteSizeHistory<NumberType> & g_rec,
536  const std::pair<NumberType, NumberType> bounds)> &choose,
537  const NumberType a1,
538  const NumberType eta,
539  const NumberType mu,
540  const NumberType a_max,
541  const unsigned int max_evaluations,
542  const bool debug_output)
543  {
544  // Note that scipy use dcsrch() from Minpack2 Fortran lib for line search
545  Assert(mu < 0.5 && mu > 0, ExcMessage("mu is not in (0,1/2)."));
546  Assert(eta < 1. && eta > mu, ExcMessage("eta is not in (mu,1)."));
547  Assert(a_max > 0, ExcMessage("max is not positive."));
548  Assert(a1 > 0 && a1 <= a_max, ExcMessage("a1 is not in (0,max]."));
549  Assert(g0 < 0, ExcMessage("Initial slope is not negative"));
550 
551  // Growth parameter for bracketing phase:
552  // 1 < tau1
553  const NumberType tau1 = 9.;
554  // shrink parameters for sectioning phase to prevent ai from being
555  // arbitrary close to the extremes of the interval.
556  // 0 < tau2 < tau3 <= 1/2
557  // tau2 <= eta is advisable
558  const NumberType tau2 = 0.1; // bound for closeness to a_lo
559  const NumberType tau3 = 0.5; // bound for closeness to a_hi
560 
561  const NumberType g0_abs = std::abs(g0);
562  const NumberType f_min = f0 + a_max * mu * g0;
563 
564  // return True if the first Wolfe condition (sufficient decrease) is
565  // satisfied
566  const auto w1 = [&](const NumberType a, const NumberType f) {
567  return f <= f0 + a * mu * g0;
568  };
569 
570  // return True if the second Wolfe condition (curvature condition) is
571  // satisfied
572  const auto w2 = [&](const NumberType g) {
573  return std::abs(g) <= eta * g0_abs;
574  };
575 
576  // Bracketing phase (Algorithm 2.6.2): look for a non-trivial interval
577  // which is known to contain an interval of acceptable points.
578  // We adopt notation of Noceal.
579  const NumberType x = std::numeric_limits<NumberType>::signaling_NaN();
580  NumberType a_lo = x, f_lo = x, g_lo = x;
581  NumberType a_hi = x, f_hi = x, g_hi = x;
582  NumberType ai = x, fi = x, gi = x;
583 
584  // count function calls in i:
585  unsigned int i = 0;
586  {
587  NumberType f_prev, g_prev, a_prev;
588  ai = a1;
589  f_prev = f0;
590  g_prev = g0;
591  a_prev = 0;
592 
593  while (i < max_evaluations)
594  {
595  const auto fgi = func(ai);
596  fi = fgi.first;
597  gi = fgi.second;
598  i++;
599 
600  if (debug_output)
601  deallog << "Bracketing phase: " << i << std::endl
602  << ai << " " << fi << " " << gi << " " << w1(ai, fi) << " "
603  << w2(gi) << " " << f_min << std::endl;
604 
605  // first check if we can stop bracketing or the whole line search:
606  if (fi <= f_min || ai == a_max)
607  {
608  if (debug_output)
609  deallog << "Reached the maximum step size." << std::endl;
610  return std::make_pair(ai, i);
611  }
612 
613  if (!w1(ai, fi) ||
614  (fi >= f_prev && i > 1)) // violate first Wolfe or not descending
615  {
616  a_lo = a_prev;
617  f_lo = f_prev;
618  g_lo = g_prev;
619 
620  a_hi = ai;
621  f_hi = fi;
622  g_hi = gi;
623  break; // end bracketing
624  }
625 
626  if (w2(gi)) // satisfies both Wolfe conditions
627  {
628  if (debug_output)
629  deallog << "Satisfied both Wolfe conditions during Bracketing."
630  << std::endl;
631 
632  Assert(w1(ai, fi), ExcInternalError());
633  return std::make_pair(ai, i);
634  }
635 
636  if (gi >= 0) // not descending
637  {
638  a_lo = ai;
639  f_lo = fi;
640  g_lo = gi;
641 
642  a_hi = a_prev;
643  f_hi = f_prev;
644  g_hi = g_prev;
645  break; // end bracketing
646  }
647 
648  // extrapolation step with the bounds
649  const auto bounds =
650  std::make_pair(2. * ai - a_prev,
651  std::min(a_max, ai + tau1 * (ai - a_prev)));
652 
653  a_prev = ai;
654  f_prev = fi;
655  g_prev = gi;
656 
657  // NOTE: Fletcher's 2.6.2 includes optional extrapolation, we
658  // simply take the upper bound
659  // Scipy increases by factor of two:
660  // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L447
661  ai = bounds.second;
662  }
663  }
664 
665  AssertThrow(
666  i < max_evaluations,
667  ExcMessage(
668  "Could not find the initial bracket within the given number of iterations."));
669 
670  // Check properties of the bracket (Theorem 3.2 in More and Thuente, 94
671  // and Eq. 2.6.3 in Fletcher 2013
672 
673  // FIXME: these conditions are actually violated for Fig3 and a1=10^3 in
674  // More and Thorenton, 94.
675 
676  /*
677  Assert((f_lo < f_hi) && w1(a_lo, f_lo), ExcInternalError());
678  Assert(((a_hi - a_lo) * g_lo < 0) && !w2(g_lo), ExcInternalError());
679  Assert((w1(a_hi, f_hi) || f_hi >= f_lo), ExcInternalError());
680  */
681 
682  // keep short history of last points to improve interpolation
683  FiniteSizeHistory<NumberType> a_rec(5), f_rec(5), g_rec(5);
684  // if neither a_lo nor a_hi are zero:
685  if (std::abs(a_lo) > std::numeric_limits<NumberType>::epsilon() &&
686  std::abs(a_hi) > std::numeric_limits<NumberType>::epsilon())
687  {
688  a_rec.add(0);
689  f_rec.add(f0);
690  g_rec.add(g0);
691  }
692 
693  // Now sectioning phase: we allow both [a_lo, a_hi] and [a_hi, a_lo]
694  while (i < max_evaluations)
695  {
696  const NumberType a_lo_safe = a_lo + tau2 * (a_hi - a_lo);
697  const NumberType a_hi_safe = a_hi - tau3 * (a_hi - a_lo);
698  const auto bounds = std::minmax(a_lo_safe, a_hi_safe);
699 
700  ai = choose(
701  a_lo, f_lo, g_lo, a_hi, f_hi, g_hi, a_rec, f_rec, g_rec, bounds);
702 
703  const std::pair<NumberType, NumberType> fgi = func(ai);
704  fi = fgi.first;
705  gi = fgi.second;
706  i++;
707 
708  if (debug_output)
709  deallog << "Sectioning phase: " << i << std::endl
710  << a_lo << " " << f_lo << " " << g_lo << " " << w1(a_lo, f_lo)
711  << " " << w2(g_lo) << std::endl
712  << a_hi << " " << f_hi << " " << g_hi << " " << w1(a_hi, f_hi)
713  << " " << w2(g_hi) << std::endl
714  << ai << " " << fi << " " << gi << " " << w1(ai, fi) << " "
715  << w2(gi) << std::endl;
716 
717  if (!w1(ai, fi) || fi >= f_lo)
718  // take [a_lo, ai]
719  {
720  a_rec.add(a_hi);
721  f_rec.add(f_hi);
722  g_rec.add(g_hi);
723 
724  a_hi = ai;
725  f_hi = fi;
726  g_hi = gi;
727  }
728  else
729  {
730  if (w2(gi)) // satisfies both wolf
731  {
732  if (debug_output)
733  deallog << "Satisfied both Wolfe conditions." << std::endl;
734  Assert(w1(ai, fi), ExcInternalError());
735  return std::make_pair(ai, i);
736  }
737 
738  if (gi * (a_hi - a_lo) >= 0)
739  // take [ai, a_lo]
740  {
741  a_rec.add(a_hi);
742  f_rec.add(f_hi);
743  g_rec.add(g_hi);
744 
745  a_hi = a_lo;
746  f_hi = f_lo;
747  g_hi = g_lo;
748  }
749  else
750  // take [ai, a_hi]
751  {
752  a_rec.add(a_lo);
753  f_rec.add(f_lo);
754  g_rec.add(g_lo);
755  }
756 
757  a_lo = ai;
758  f_lo = fi;
759  g_lo = gi;
760  }
761  }
762 
763  // if we got here, we could not find the solution
764  AssertThrow(
765  false,
766  ExcMessage(
767  "Could not could complete the sectioning phase within the given number of iterations."));
768  return std::make_pair(std::numeric_limits<NumberType>::signaling_NaN(), i);
769  }
770 
771 #endif // DOXYGEN
772 
773 } // namespace LineMinimization
774 
775 DEAL_II_NAMESPACE_CLOSE
776 
777 #endif // dealii_line_minimization_h
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::size_t size() const
boost::optional< NumberType > cubic_fit_three_points(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType x_rec, const NumberType f_rec)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
boost::optional< NumberType > cubic_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi)
static ::ExceptionBase & ExcMessage(std::string arg1)
NumberType poly_fit_three_points(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)
#define Assert(cond, exc)
Definition: exceptions.h:1407
NumberType poly_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)
void add(const T &element)
boost::optional< NumberType > quadratic_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi)
std::pair< NumberType, unsigned int > line_search(const std::function< std::pair< NumberType, NumberType >(const NumberType x)> &func, const NumberType f0, const NumberType g0, const std::function< NumberType(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)> &interpolate, const NumberType a1, const NumberType eta=0.9, const NumberType mu=0.01, const NumberType a_max=std::numeric_limits< NumberType >::max(), const unsigned int max_evaluations=20, const bool debug_output=false)
static ::ExceptionBase & ExcInternalError()