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