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