Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_line_minimization_h
16#define dealii_line_minimization_h
17
18#include <deal.II/base/config.h>
19
24
26
27#include <algorithm>
28#include <fstream>
29#include <functional>
30#include <limits>
31#include <optional>
32#include <string>
33
34
36
41{
50 template <typename NumberType>
51 std::optional<NumberType>
52 quadratic_fit(const NumberType x_low,
53 const NumberType f_low,
54 const NumberType g_low,
55 const NumberType x_hi,
56 const NumberType f_hi);
57
67 template <typename NumberType>
68 std::optional<NumberType>
69 cubic_fit(const NumberType x_low,
70 const NumberType f_low,
71 const NumberType g_low,
72 const NumberType x_hi,
73 const NumberType f_hi,
74 const NumberType g_hi);
75
83 template <typename NumberType>
84 std::optional<NumberType>
85 cubic_fit_three_points(const NumberType x_low,
86 const NumberType f_low,
87 const NumberType g_low,
88 const NumberType x_hi,
89 const NumberType f_hi,
90 const NumberType x_rec,
91 const NumberType f_rec);
92
104 template <typename NumberType>
105 NumberType
106 poly_fit(const NumberType x_low,
107 const NumberType f_low,
108 const NumberType g_low,
109 const NumberType x_hi,
110 const NumberType f_hi,
111 const NumberType g_hi,
115 const std::pair<NumberType, NumberType> bounds);
116
121 template <typename NumberType>
122 NumberType
123 poly_fit_three_points(const NumberType x_low,
124 const NumberType f_low,
125 const NumberType g_low,
126 const NumberType x_hi,
127 const NumberType f_hi,
128 const NumberType g_hi,
132 const std::pair<NumberType, NumberType> bounds);
133
134
321 template <typename NumberType>
322 std::pair<NumberType, unsigned int>
324 const std::function<std::pair<NumberType, NumberType>(const NumberType x)>
325 &func,
326 const NumberType f0,
327 const NumberType g0,
328 const std::function<
329 NumberType(const NumberType x_low,
330 const NumberType f_low,
331 const NumberType g_low,
332 const NumberType x_hi,
333 const NumberType f_hi,
334 const NumberType g_hi,
338 const std::pair<NumberType, NumberType> bounds)> &interpolate,
339 const NumberType a1,
340 const NumberType eta = 0.9,
341 const NumberType mu = 0.01,
342 const NumberType a_max = std::numeric_limits<NumberType>::max(),
343 const unsigned int max_evaluations = 20,
344 const bool debug_output = false);
345
346
347 // ------------------- inline and template functions ----------------
348
349
350#ifndef DOXYGEN
351
352
353 template <typename NumberType>
354 std::optional<NumberType>
355 quadratic_fit(const NumberType x1,
356 const NumberType f1,
357 const NumberType g1,
358 const NumberType x2,
359 const NumberType f2)
360 {
361 Assert(x1 != x2, ExcMessage("Point are the same"));
362 const NumberType denom = (2. * g1 * x2 - 2. * g1 * x1 - 2. * f2 + 2. * f1);
363 if (denom == 0)
364 return {};
365 else
366 return (g1 * (x2 * x2 - x1 * x1) + 2. * (f1 - f2) * x1) / denom;
367 }
368
369
370
371 template <typename NumberType>
372 std::optional<NumberType>
373 cubic_fit(const NumberType x1,
374 const NumberType f1,
375 const NumberType g1,
376 const NumberType x2,
377 const NumberType f2,
378 const NumberType g2)
379 {
380 Assert(x1 != x2, ExcMessage("Points are the same"));
381 const NumberType beta1 = g1 + g2 - 3. * (f1 - f2) / (x1 - x2);
382 const NumberType s = beta1 * beta1 - g1 * g2;
383 if (s < 0)
384 return {};
385
386 const NumberType beta2 = std::sqrt(s);
387 const NumberType denom =
388 x1 < x2 ? g2 - g1 + 2. * beta2 : g1 - g2 + 2. * beta2;
389 if (denom == 0.)
390 return {};
391
392 return x1 < x2 ? x2 - (x2 - x1) * (g2 + beta2 - beta1) / denom :
393 x1 - (x1 - x2) * (g1 + beta2 - beta1) / denom;
394 }
395
396
397
398 template <typename NumberType>
399 std::optional<NumberType>
400 cubic_fit_three_points(const NumberType x1,
401 const NumberType f1,
402 const NumberType g1,
403 const NumberType x2,
404 const NumberType f2,
405 const NumberType x3,
406 const NumberType f3)
407 {
408 Assert(x1 != x2, ExcMessage("Points are the same"));
409 Assert(x1 != x3, ExcMessage("Points are the same"));
410 // f(x) = A *(x-x1)^3 + B*(x-x1)^2 + C*(x-x1) + D
411 // =>
412 // D = f1
413 // C = g1
414
415 // the rest is a system of 2 equations:
416
417 const NumberType x2_shift = x2 - x1;
418 const NumberType x3_shift = x3 - x1;
419 const NumberType r1 = f2 - f1 - g1 * x2_shift;
420 const NumberType r2 = f3 - f1 - g1 * x3_shift;
421 const NumberType denom =
422 Utilities::fixed_power<2>(x2_shift * x3_shift) * (x2_shift - x3_shift);
423 if (denom == 0.)
424 return {};
425
426 const NumberType A = (r1 * Utilities::fixed_power<2>(x3_shift) -
427 r2 * Utilities::fixed_power<2>(x2_shift)) /
428 denom;
429 const NumberType B = (r2 * Utilities::fixed_power<3>(x2_shift) -
430 r1 * Utilities::fixed_power<3>(x3_shift)) /
431 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 {};
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 std::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*/,
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 std::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 std::optional<NumberType>{};
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>
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,
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
666 i < max_evaluations,
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
765 false,
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
776
777#endif // dealii_line_minimization_h
std::size_t size() const
void add(const T &element)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:36
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)
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::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::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)
std::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)
STL namespace.
::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 > &)