From: Denis Davydov Date: Thu, 13 Sep 2018 10:14:53 +0000 (+0200) Subject: add line minimization functions X-Git-Tag: v9.1.0-rc1~696^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=26d62587244b9e39f124f587471d87c692f34634;p=dealii.git add line minimization functions --- diff --git a/doc/news/changes/minor/20180913DenisDavydov b/doc/news/changes/minor/20180913DenisDavydov new file mode 100644 index 0000000000..31324ece97 --- /dev/null +++ b/doc/news/changes/minor/20180913DenisDavydov @@ -0,0 +1,3 @@ +New: Add line minimization functions. +
+(Denis Davydov 2018/09/13) diff --git a/include/deal.II/optimization/line_minimization.h b/include/deal.II/optimization/line_minimization.h new file mode 100644 index 0000000000..160913b74d --- /dev/null +++ b/include/deal.II/optimization/line_minimization.h @@ -0,0 +1,606 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +#ifndef dealii_line_minimization_h +#define dealii_line_minimization_h + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +using namespace dealii; + +/** + * A namespace for various algorithms related to minimization over line. + */ +namespace LineMinimization +{ + /** + * Given $x\_low$ and $x\_hi$ together with values of function + * $f(x\_low)$ and $f(x\_hi)$ and the gradient $g(x\_low)$, return the local + * minimizer of the quadratic interpolation function. + * + * The return type is optional to fit with similar function which may + * not have solution for given parameters. + */ + template + boost::optional + quadratic_fit(const NumberType x_low, + const NumberType f_low, + const NumberType g_low, + const NumberType x_hi, + const NumberType f_hi); + + /** + * Given $x\_low$ and $x\_hi$ together with values of function + * $f(x\_low)$ and $f(x\_hi)$) and its gradients ($g(x\_low)*g(x\_hi) < 0$) at + * those points, return the local minimizer of the cubic interpolation + * function. That is the location where the cubic interpolation function + * attains its minimum. + * + * The return type is optional as the real-valued solution might not exist. + */ + template + boost::optional + 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); + + /** + * Find the minimizer of a cubic polynomial that goes through the + * points $f\_low=f(x\_low)$, $f\_hi=f(x\_hi)$ and $f\_rec(x\_rec)$ + * and has derivatve $g\_low$ at $x\_low$. + * + * The return type is optional as the real-valued solution might not exist. + */ + template + boost::optional + 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); + + /** + * Return the minimizer of a polynomial using function values @p f_low @p f_hi @p f_rec[0] at three points + * @p x_low @p x_hi @p x_rec[0] and derivatives at two points @p g_low and @p g_hi. The returned point + * should be within the bounds @p bounds . + * + * This function will first try the cubic fit (see cubic_fit ). If it's + * unsuccessfull or not + * within the provided @p bounds the quadratic fit will be performed (see quadratic_fit ). The function will + * fallback to bisection if quadratic fit fails. + */ + template + 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 & x_rec, + const FiniteSizeHistory & f_rec, + const FiniteSizeHistory & g_rec, + const std::pair bounds); + + /** + * Same as above but doing cubic fit with three points (see + * cubic_fit_three_points ). + */ + template + 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 & x_rec, + const FiniteSizeHistory & f_rec, + const FiniteSizeHistory & g_rec, + const std::pair bounds); + + + /** + * Perform a line search in $(0,max]$ with strong Wolfe conditions + * \f[ + * f(\alpha) \le f(0) + \alpha \mu f'(0) \\ + * |f'(\alpha)| \le \eta |f'(0)| + * \f] + * using one dimensional + * functions @p func and a function @p interpolate to choose a new point + * from the interval based on the function values and derivatives at its ends. + * @p a1 is a trial estimate of the first step. + * Interpolation can be done using poly_fit or poly_fit_three_points . + * + * The function implements Algorithms + * 2.6.2 and 2.6.4 on pages 34-35 in Fletcher, 2013, Practical methods of + * optimization. these are minor variations of the Algorithm 3.5 and 3.6 on + * pages 60-61 in Nocedal and Wright, Numerical optimization. + * It consists of a bracketing phase and a zoom phase, where @p interpolate is used. + * + * The function returns the step size and the number of times function @p func was called. + * + * @param func a one dimensional function which returns value and derivative + * at the given point. + * @param f0 function value the origin + * @param g0 function derivative the origin + * @param interpolate a function which determines how interpolation is done + * during the zoom phase. It takes values and derivatives at the current + * interval/bracket ($f\_low$, $f\_hi$) as well as up to 5 values and + * derivatives at previous steps. The returned value is to be provided within + * the given bounds. + * @param a1 initial trial step for bracketing phase + * @param eta a parameter in the second Wolfe condition (curvature condition) + * @param mu a parameter in the first Wolfe condition (sufficient decrease) + * @param a_max maximum allowed step size + * @param max_evaluations maximum allowed number of function evaluations + * @param debug_output a flag to do extra debug output into deallog static + * object + */ + template + std::pair + line_search( + const std::function(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 & x_rec, + const FiniteSizeHistory & f_rec, + const FiniteSizeHistory & g_rec, + const std::pair bounds)> &interpolate, + const NumberType a1, + const NumberType eta = 0.9, + const NumberType mu = 0.01, + const NumberType a_max = std::numeric_limits::max(), + const unsigned int max_evaluations = 20, + const bool debug_output = false); + + // ------------------- inline and template functions ---------------- +#ifndef DOXYGEN + + template + boost::optional + quadratic_fit(const NumberType x1, + const NumberType f1, + const NumberType g1, + const NumberType x2, + const NumberType f2) + { + Assert(x1 != x2, ExcMessage("Point are the same")); + const NumberType denom = (2. * g1 * x2 - 2. * g1 * x1 - 2. * f2 + 2. * f1); + if (denom == 0) + return boost::none; + else + return (g1 * (x2 * x2 - x1 * x1) + 2. * (f1 - f2) * x1) / denom; + } + + template + boost::optional + cubic_fit(const NumberType x1, + const NumberType f1, + const NumberType g1, + const NumberType x2, + const NumberType f2, + const NumberType g2) + { + Assert(x1 != x2, ExcMessage("Points are the same")); + const NumberType beta1 = g1 + g2 - 3. * (f1 - f2) / (x1 - x2); + const NumberType s = beta1 * beta1 - g1 * g2; + if (s < 0) + return boost::none; + + const NumberType beta2 = std::sqrt(s); + const NumberType denom = + x1 < x2 ? g2 - g1 + 2. * beta2 : g1 - g2 + 2. * beta2; + if (denom == 0.) + return boost::none; + + return x1 < x2 ? x2 - (x2 - x1) * (g2 + beta2 - beta1) / denom : + x1 - (x1 - x2) * (g1 + beta2 - beta1) / denom; + } + + + + template + boost::optional + cubic_fit_three_points(const NumberType x1, + const NumberType f1, + const NumberType g1, + const NumberType x2, + const NumberType f2, + const NumberType x3, + const NumberType f3) + { + Assert(x1 != x2, ExcMessage("Points are the same")); + Assert(x1 != x3, ExcMessage("Points are the same")); + // f(x) = A *(x-x1)^3 + B*(x-x1)^2 + C*(x-x1) + D + // => + // D = f1 + // C = g1 + + // the rest is a system of 2 equations: + + const NumberType x2_shift = x2 - x1; + const NumberType x3_shift = x3 - x1; + const NumberType r1 = f2 - f1 - g1 * x2_shift; + const NumberType r2 = f3 - f1 - g1 * x3_shift; + const NumberType denom = + std::pow(x2_shift * x3_shift, 2) * (x2_shift - x3_shift); + if (denom == 0.) + return boost::none; + + const NumberType A = + (r1 * std::pow(x3_shift, 2) - r2 * std::pow(x2_shift, 2)) / denom; + const NumberType B = + (r2 * std::pow(x2_shift, 3) - r1 * std::pow(x3_shift, 3)) / denom; + const NumberType &C = g1; + + // now get the minimizer: + const NumberType radical = B * B - A * C * 3; + if (radical < 0) + return boost::none; + + return x1 + (-B + std::sqrt(radical)) / (A * 3); + } + + + template + NumberType + poly_fit(const NumberType x1, + const NumberType f1, + const NumberType g1, + const NumberType x2, + const NumberType f2, + const NumberType g2, + const FiniteSizeHistory &, + const FiniteSizeHistory &, + const FiniteSizeHistory &, + const std::pair bounds) + { + Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds")); + + // Similar to scipy implementation but we fit based on two points + // with their gradients and do bisection on bounds. + // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563 + + // First try cubic interpolation + boost::optional res = cubic_fit(x1, f1, g1, x2, f2, g2); + if (res && *res >= bounds.first && *res <= bounds.second) + return *res; + + // cubic either fails or outside of safe region, do quadratic: + res = quadratic_fit(x1, f1, g1, x2, f2); + if (res && *res >= bounds.first && *res <= bounds.second) + return *res; + + // quadratic either failed or outside of safe region. Do bisection + // on safe region + return (bounds.first + bounds.second) * 0.5; + } + + + + template + NumberType + poly_fit_three_points(const NumberType x1, + const NumberType f1, + const NumberType g1, + const NumberType x2, + const NumberType f2, + const NumberType g2, + const FiniteSizeHistory &x_rec, + const FiniteSizeHistory &f_rec, + const FiniteSizeHistory & /*g_rec*/, + const std::pair bounds) + { + Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds")); + AssertDimension(x_rec.size(), f_rec.size()); + + // Same as scipy implementation where cubic fit is using 3 points + // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563 + + // First try cubic interpolation after first iteration + boost::optional res = + x_rec.size() > 0 ? + cubic_fit_three_points(x1, f1, g1, x2, f2, x_rec[0], f_rec[0]) : + boost::none; + if (res && *res >= bounds.first && *res <= bounds.second) + return *res; + + // cubic either fails or outside of safe region, do quadratic: + res = quadratic_fit(x1, f1, g1, x2, f2); + if (res && *res >= bounds.first && *res <= bounds.second) + return *res; + + // quadratic either failed or outside of safe region. Do bisection + // on safe region + return (bounds.first + bounds.second) * 0.5; + } + + + template + std::pair + line_search( + const std::function(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 & x_rec, + const FiniteSizeHistory & f_rec, + const FiniteSizeHistory & g_rec, + const std::pair bounds)> &choose, + const NumberType a1, + const NumberType eta, + const NumberType mu, + const NumberType a_max, + const unsigned int max_evaluations, + const bool debug_output) + { + // Note that scipy use dcsrch() from Minpack2 Fortran lib for line search + Assert(mu < 0.5 && mu > 0, ExcMessage("mu is not in (0,1/2).")); + Assert(eta < 1. && eta > mu, ExcMessage("eta is not in (mu,1).")); + Assert(a_max > 0, ExcMessage("max is not positive.")); + Assert(a1 > 0 && a1 <= a_max, ExcMessage("a1 is not in (0,max].")); + Assert(g0 < 0, ExcMessage("Initial slope is not negative")); + + // Growth parameter for bracketing phase: + // 1 < tau1 + const NumberType tau1 = 9.; + // shrink parameters for sectioning phase to prevent ai from being + // arbitrary close to the extremes of the interval. + // 0 < tau2 < tau3 <= 1/2 + // tau2 <= eta is advisable + const NumberType tau2 = 0.1; // bound for closeness to a_lo + const NumberType tau3 = 0.5; // bound for closeness to a_hi + + const NumberType g0_abs = std::abs(g0); + const NumberType f_min = f0 + a_max * mu * g0; + + // return True if the first Wolfe condition (sufficient decrease) is + // satisfied + const auto w1 = [&](const NumberType a, const NumberType f) { + return f <= f0 + a * mu * g0; + }; + + // return True if the second Wolfe condition (curvature condition) is + // satisfied + const auto w2 = [&](const NumberType g) { + return std::abs(g) <= eta * g0_abs; + }; + + // Bracketing phase (Algorithm 2.6.2): look for a non-trivial interval + // which is known to contain an interval of acceptable points. + // We adopt notation of Noceal. + const NumberType x = std::numeric_limits::signaling_NaN(); + NumberType a_lo = x, f_lo = x, g_lo = x; + NumberType a_hi = x, f_hi = x, g_hi = x; + NumberType ai = x, fi = x, gi = x; + + // count function calls in i: + unsigned int i = 0; + { + NumberType f_prev, g_prev, a_prev; + ai = a1; + f_prev = f0; + g_prev = g0; + a_prev = 0; + + while (i < max_evaluations) + { + const auto fgi = func(ai); + fi = fgi.first; + gi = fgi.second; + i++; + + if (debug_output) + deallog << "Bracketing phase: " << i << std::endl + << ai << " " << fi << " " << gi << std::endl; + + // first check if we can stop bracketing or the whole line search: + if (fi <= f_min || ai == a_max) + return std::make_pair(ai, i); + + if (!w1(ai, fi) || + (fi >= f_prev && i > 1)) // violate first Wolfe or not descending + { + a_lo = a_prev; + f_lo = f_prev; + g_lo = g_prev; + + a_hi = ai; + f_hi = fi; + g_hi = gi; + break; // end bracketing + } + + if (w2(gi)) // satisfies both Wolfe conditions + { + Assert(w1(ai, fi), ExcInternalError()); + return std::make_pair(ai, i); + } + + if (gi >= 0) // not descending + { + a_lo = ai; + f_lo = fi; + g_lo = gi; + + a_hi = a_prev; + f_hi = f_prev; + g_hi = g_prev; + break; // end bracketing + } + + // extrapolation step with the bounds + const auto bounds = + std::make_pair(2. * ai - a_prev, + std::min(a_max, ai + tau1 * (ai - a_prev))); + + a_prev = ai; + f_prev = fi; + g_prev = gi; + + // NOTE: Fletcher's 2.6.2 includes optional extrapolation, we + // simply take the upper bound + // Scipy increases by factor of two: + // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L447 + ai = bounds.second; + } + } + + AssertThrow( + i < max_evaluations, + ExcMessage( + "Could not find the initial bracket within the given number of iterations.")); + + // Check properties of the bracket (Theorem 3.2 in More and Thuente, 94 + // and Eq. 2.6.3 in Fletcher 2013 + + // FIXME: these conditions are actually violated for Fig3 and a1=10^3 in + // More and Thorenton, 94. + + /* + Assert((f_lo < f_hi) && w1(a_lo, f_lo), ExcInternalError()); + Assert(((a_hi - a_lo) * g_lo < 0) && !w2(g_lo), ExcInternalError()); + Assert((w1(a_hi, f_hi) || f_hi >= f_lo), ExcInternalError()); + */ + + // keep short history of last points to improve interpolation + FiniteSizeHistory a_rec(5), f_rec(5), g_rec(5); + // if neither a_lo nor a_hi are zero: + if (std::abs(a_lo) > std::numeric_limits::epsilon() && + std::abs(a_hi) > std::numeric_limits::epsilon()) + { + a_rec.add(0); + f_rec.add(f0); + g_rec.add(g0); + } + + // Now sectioning phase: we allow both [a_lo, a_hi] and [a_hi, a_lo] + while (i < max_evaluations) + { + const NumberType a_lo_safe = a_lo + tau2 * (a_hi - a_lo); + const NumberType a_hi_safe = a_hi - tau3 * (a_hi - a_lo); + const auto bounds = std::minmax(a_lo_safe, a_hi_safe); + + ai = choose( + a_lo, f_lo, g_lo, a_hi, f_hi, g_hi, a_rec, f_rec, g_rec, bounds); + + const std::pair fgi = func(ai); + fi = fgi.first; + gi = fgi.second; + i++; + + if (debug_output) + deallog << "Sectioning phase: " << i << std::endl + << a_lo << " " << f_lo << " " << g_lo << std::endl + << a_hi << " " << f_hi << " " << g_hi << std::endl + << ai << " " << fi << " " << gi << std::endl; + + if (!w1(ai, fi) || fi >= f_lo) + // take [a_lo, ai] + { + a_rec.add(a_hi); + f_rec.add(f_hi); + g_rec.add(g_hi); + + a_hi = ai; + f_hi = fi; + g_hi = gi; + } + else + { + if (w2(gi)) // satisfies both wolf + { + Assert(w1(ai, fi), ExcInternalError()); + return std::make_pair(ai, i); + } + + if (gi * (a_hi - a_lo) >= 0) + // take [ai, a_lo] + { + a_rec.add(a_hi); + f_rec.add(f_hi); + g_rec.add(g_hi); + + a_hi = a_lo; + f_hi = f_lo; + g_hi = g_lo; + } + else + // take [ai, a_hi] + { + a_rec.add(a_lo); + f_rec.add(f_lo); + g_rec.add(g_lo); + } + + a_lo = ai; + f_lo = fi; + g_lo = gi; + } + } + + // if we got here, we could not find the solution + AssertThrow( + false, + ExcMessage( + "Could not could complete the sectioning phase within the given number of iterations.")); + return std::make_pair(std::numeric_limits::signaling_NaN(), i); + } + +#endif + +} // namespace LineMinimization + +DEAL_II_NAMESPACE_CLOSE + +#endif // dealii_line_minimization_h diff --git a/tests/optimization/CMakeLists.txt b/tests/optimization/CMakeLists.txt new file mode 100644 index 0000000000..3772facf43 --- /dev/null +++ b/tests/optimization/CMakeLists.txt @@ -0,0 +1,5 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) +INCLUDE(../setup_testsubproject.cmake) +PROJECT(testsuite CXX) +DEAL_II_PICKUP_TESTS() + diff --git a/tests/optimization/cubic_fit.cc b/tests/optimization/cubic_fit.cc new file mode 100644 index 0000000000..742fa702d6 --- /dev/null +++ b/tests/optimization/cubic_fit.cc @@ -0,0 +1,67 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check minimization of the cubic fit based on f(x1), f(x2) and f'(x1) +// and f'(x2) + +#include + +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +void +test() +{ + // test 1: + { + auto f = [](double x) { + return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x; + }; + auto g = [](double x) { + return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1; + }; + + const double x1 = 5; + const double x2 = 17; + const double f1 = f(x1); + const double f2 = f(x2); + const double g1 = g(x1); + const double g2 = g(x2); + const double res = *LineMinimization::cubic_fit(x1, f1, g1, x2, f2, g2); + const double res2 = *LineMinimization::cubic_fit(x2, f2, g2, x1, f1, g1); + deallog << x1 << " " << x2 << std::endl + << f1 << " " << f2 << std::endl + << g1 << " " << g2 << std::endl + << res << std::endl + << res2 << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/cubic_fit.output b/tests/optimization/cubic_fit.output new file mode 100644 index 0000000000..b0ad6960d9 --- /dev/null +++ b/tests/optimization/cubic_fit.output @@ -0,0 +1,5 @@ +DEAL::5.00000 17.0000 +DEAL::-1874.50 -14737.3 +DEAL::-999.900 2312.10 +DEAL::14.6115 +DEAL::14.6115 diff --git a/tests/optimization/cubic_fit_three_points.cc b/tests/optimization/cubic_fit_three_points.cc new file mode 100644 index 0000000000..0f970ec6c1 --- /dev/null +++ b/tests/optimization/cubic_fit_three_points.cc @@ -0,0 +1,67 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check minimization of the cubic fit based on f(x1), f(x2) and f'(x1) +// and f'(x2) + +#include + +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +void +test() +{ + // test 1: + { + auto f = [](double x) { + return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x; + }; + auto g = [](double x) { + return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1; + }; + + const double x1 = 17; + const double x2 = 10; + const double x3 = 5; + const double f1 = f(x1); + const double f2 = f(x2); + const double f3 = f(x3); + const double g1 = g(x1); + const double res = + *LineMinimization::cubic_fit_three_points(x1, f1, g1, x2, f2, x3, f3); + deallog << x1 << " " << f1 << " " << g1 << std::endl + << x2 << " " << f2 << std::endl + << x3 << " " << f3 << std::endl + << res << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/cubic_fit_three_points.output b/tests/optimization/cubic_fit_three_points.output new file mode 100644 index 0000000000..fa97885b14 --- /dev/null +++ b/tests/optimization/cubic_fit_three_points.output @@ -0,0 +1,4 @@ +DEAL::17.0000 -14737.3 2312.10 +DEAL::10.0000 -9999.00 +DEAL::5.00000 -1874.50 +DEAL::14.8441 diff --git a/tests/optimization/line_minimization.cc b/tests/optimization/line_minimization.cc new file mode 100644 index 0000000000..8dfb886e62 --- /dev/null +++ b/tests/optimization/line_minimization.cc @@ -0,0 +1,83 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check line minimization with strong Wolfe conditions + +#include + +#include + +#include +#include + +#include "../tests.h" +using namespace dealii; + +/* + * MWE in Maxima + +func(x):=100 * x^4 + (1-x)^2; +gfunc(x):=''(diff(func(x),x)); +eta : 0.1; +mu : 0.01; +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +plot2d([func(x),gfunc(x),w1(x), w2(x)], [x,0.1,0.2]); +w1(0.159668); +w2(0.159668); +bfloat(solve(gfunc(x)=0)[3]); + + */ + + +void +test() +{ + // test 1: + { + const double min_x = 0.161262023139589; + auto func = [](const double x) { + const double f = 100. * std::pow(x, 4) + std::pow(1. - x, 2); + const double g = 400. * std::pow(x, 3) - 2. * (1. - x); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + const auto res = + LineMinimization::line_search(func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + 0.1, + 0.1, + 0.01, + 100, + 20, + true); + deallog << "Solution: " << res.first << std::endl + << "Distance: " << std::fabs(res.first - min_x) << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/line_minimization.output b/tests/optimization/line_minimization.output new file mode 100644 index 0000000000..e1a7ffbcd4 --- /dev/null +++ b/tests/optimization/line_minimization.output @@ -0,0 +1,18 @@ +DEAL::Bracketing phase: 1 +DEAL::0.100000 0.820000 -1.40000 +DEAL::Bracketing phase: 2 +DEAL::1.00000 100.000 400.000 +DEAL::Sectioning phase: 3 +DEAL::0.100000 0.820000 -1.40000 +DEAL::1.00000 100.000 400.000 +DEAL::0.333333 1.67901 13.4815 +DEAL::Sectioning phase: 4 +DEAL::0.100000 0.820000 -1.40000 +DEAL::0.333333 1.67901 13.4815 +DEAL::0.167641 0.771802 0.219785 +DEAL::Sectioning phase: 5 +DEAL::0.167641 0.771802 0.219785 +DEAL::0.100000 0.820000 -1.40000 +DEAL::0.159668 0.771152 -0.0524435 +DEAL::Solution: 0.159668 +DEAL::Distance: 0.00159407 diff --git a/tests/optimization/line_minimization_02.cc b/tests/optimization/line_minimization_02.cc new file mode 100644 index 0000000000..cc49a4f97a --- /dev/null +++ b/tests/optimization/line_minimization_02.cc @@ -0,0 +1,133 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check line minimization with strong Wolfe conditions +// similar to line_minimization.cc but a different function and +// different initial steps + +#include + +#include + +#include +#include + +#include "../tests.h" +using namespace dealii; + +/* + * MWE in Maxima + +func(x):=-3*x/(x^2+2)-0.03*x; +gfunc(x):=''(diff(func(x),x)); +mu:0.025; +eta : 0.1; +a_max : 30; +f_min : func(0) + a_max * mu * gfunc(0); +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 ); +pline(x):=func(0)+x*mu*gfunc(0); +plot2d([func(x),pline(x),w(x)], [x,0,20]); +bfloat(solve(gfunc(x)=0)[4]); + + */ + + +void +test() +{ + // test 1: + { + const double min_x = 1.474531468108294; + auto func = [](const double x) { + const double f = (-3. * x) / (x * x + 2.) - 0.03 * x; + const double g = + -3. / (x * x + 2.) + (6. * x * x) / std::pow(x * x + 2., 2) - 0.03; + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + { + deallog << "Case 1:" << std::endl; + // First, overshoot and get to solution immediately + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + 13, + 0.1, + 0.025, + 30, + 20, + true); + deallog << "Solution: " << res.first << std::endl + << "Distance: " << std::fabs(res.first - min_x) << std::endl; + } + + { + deallog << "Case 2:" << std::endl; + // Now a small step to converge where needed: + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + 0.1, + 0.1, + 0.025, + 30, + 20, + true); + deallog << "Solution: " << res.first << std::endl + << "Distance: " << std::fabs(res.first - min_x) << std::endl; + } + + { + deallog << "Case 3:" << std::endl; + // Now do a big step so that next one in bracketing satisfies both Wolf: + // at the termination point the derivative is alos negative! + // Also that point satisfies both Wolfe conditions as well, but + // we are interested in another segment, which contains local + // minimizer + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + 1, + 0.1, + 0.025, + 30, + 20, + true); + deallog << "Solution: " << res.first << std::endl + << "Distance: " << std::fabs(res.first - min_x) << std::endl; + } + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/line_minimization_02.output b/tests/optimization/line_minimization_02.output new file mode 100644 index 0000000000..5a8e71a380 --- /dev/null +++ b/tests/optimization/line_minimization_02.output @@ -0,0 +1,37 @@ +DEAL::Case 1: +DEAL::Bracketing phase: 1 +DEAL::13.0000 -0.618070 -0.0128665 +DEAL::Solution: 13.0000 +DEAL::Distance: 11.5255 +DEAL::Case 2: +DEAL::Bracketing phase: 1 +DEAL::0.100000 -0.152254 -1.50769 +DEAL::Bracketing phase: 2 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::Bracketing phase: 3 +DEAL::9.10000 -0.594896 0.00370484 +DEAL::Sectioning phase: 4 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::9.10000 -0.594896 0.00370484 +DEAL::3.09290 -0.895024 0.139676 +DEAL::Sectioning phase: 5 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::3.09290 -0.895024 0.139676 +DEAL::1.60613 -1.10031 0.0529144 +DEAL::Solution: 1.60613 +DEAL::Distance: 0.131601 +DEAL::Case 3: +DEAL::Bracketing phase: 1 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::Bracketing phase: 2 +DEAL::10.0000 -0.594118 -0.00174164 +DEAL::Sectioning phase: 3 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::10.0000 -0.594118 -0.00174164 +DEAL::3.36365 -0.858821 0.127629 +DEAL::Sectioning phase: 4 +DEAL::1.00000 -1.03000 -0.363333 +DEAL::3.36365 -0.858821 0.127629 +DEAL::1.65166 -1.09756 0.0676998 +DEAL::Solution: 1.65166 +DEAL::Distance: 0.177131 diff --git a/tests/optimization/line_minimization_03.cc b/tests/optimization/line_minimization_03.cc new file mode 100644 index 0000000000..7774bf688d --- /dev/null +++ b/tests/optimization/line_minimization_03.cc @@ -0,0 +1,195 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check line minimization with strong Wolfe conditions +// use functions from More Thuente 1994, Line Search Algorithms with +// Guaranteed Sufficient Decrease + +#include + +#include + +#include +#include + +#include "../tests.h" +using namespace dealii; + +/* + * MWE in Maxima + + +Case 1 (compare to Table I) Function 5.1 Figure 3: + +b:2; +func(x):=-x/(x^2+b); +gfunc(x):=''(diff(func(x),x)); +mu:0.001; +eta : 0.1; +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 ); +pline(x):=func(0)+x*mu*gfunc(0); +plot2d([func(x),pline(x),w(x)], [x,0,16]); + +Case 2, Function 5.2 Figure 4 (compared to Table II we have more iterations, up +to 8 vs 16): + +b:0.004; +func(x):=(x+b)^5-2*(x+b)^4; +gfunc(x):=''(diff(func(x),x)); +mu:0.1; +eta : 0.1; +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 ); +pline(x):=func(0)+x*mu*gfunc(0); +plot2d([func(x),pline(x),w(x)], [x,0,2]); + +Case 3-5 Function 5.4 Figure 6 (Table IV-V-VI): + +b1:0.001; +b2:0.001; + +b1:0.01; +b2:0.001; + +b1:0.001; +b2:0.01; + +g(x):=(1+x^2)^(1/2)-x; +func(x):=g(b1)*((1-x)^2+b2^2)^(1/2) + g(b2)*(x^2+b1^2)^(1/2); +gfunc(x):=''(diff(func(x),x)); +plot2d([func(x)], [x,0,1]); + + */ + + +void +test() +{ + const std::vector values = {{1e-3, 1e-1, 1e+1, 1e+3}}; + double f0, g0, fi, gi; + + { + deallog << "Table 1:" << std::endl; + const double b = 2; + auto func = [&](const double &x) { + const double f = -x / (x * x + b); + const double g = 2. * x * x / std::pow(x * x + 2., 2) - 1. / (x * x + 2.); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + a1, + 0.1, + 0.001); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + + { + deallog << "Table 2:" << std::endl; + const double b = 0.004; + auto func = [&](const double x) { + const double f = std::pow(x + b, 5) - 2. * std::pow(x + b, 4); + const double g = 5. * std::pow(x + b, 4) - 8. * std::pow(x + b, 3); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + a1, + 0.100001, + 0.1); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + + { + const std::vector> params = { + {{0.001, 0.001}, {0.01, 0.001}, {0.001, 0.01}}}; + + unsigned int ind = 4; + for (auto p : params) + { + deallog << "Table " << ind++ << ":" << std::endl; + const double b1 = p.first; + const double b2 = p.second; + + const double gb1 = std::sqrt(1. + b1 * b1) - b1; + const double gb2 = std::sqrt(1. + b2 * b2) - b2; + + auto func = [&](const double x) { + const double f = gb1 * std::sqrt(std::pow(1. - x, 2) + b2 * b2) + + gb2 * std::sqrt(x * x + b1 * b1); + const double g = + gb2 * x / sqrt(x * x + b1 * b1) - + gb1 * (1. - x) / std::sqrt(std::pow(1 - x, 2) + b2 * b2); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit, + a1, + 0.00100001, + 0.001); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/line_minimization_03.output b/tests/optimization/line_minimization_03.output new file mode 100644 index 0000000000..688f24c583 --- /dev/null +++ b/tests/optimization/line_minimization_03.output @@ -0,0 +1,25 @@ +DEAL::Table 1: +DEAL::7 1.47453 0.00999934 +DEAL::5 1.51075 0.0153974 +DEAL::1 10.0000 0.00941945 +DEAL::4 37.0539 0.000725162 +DEAL::Table 2: +DEAL::19 1.59600 -2.81091e-10 +DEAL::16 1.59600 -9.10774e-11 +DEAL::13 1.59600 3.23723e-10 +DEAL::16 1.59600 2.57199e-09 +DEAL::Table 4: +DEAL::3 0.0910000 -5.97089e-05 +DEAL::1 0.100000 -4.93296e-05 +DEAL::3 0.371819 -2.34721e-06 +DEAL::5 0.545285 7.35858e-07 +DEAL::Table 5: +DEAL::5 0.0743487 3.55960e-05 +DEAL::3 0.0744514 5.98435e-05 +DEAL::10 0.0733164 -0.000213775 +DEAL::13 0.0767649 0.000581084 +DEAL::Table 6: +DEAL::12 0.925705 -2.28542e-05 +DEAL::11 0.923466 -0.000531186 +DEAL::12 0.925750 -1.23021e-05 +DEAL::9 0.924924 -0.000205208 diff --git a/tests/optimization/line_minimization_03b.cc b/tests/optimization/line_minimization_03b.cc new file mode 100644 index 0000000000..8e19d3b645 --- /dev/null +++ b/tests/optimization/line_minimization_03b.cc @@ -0,0 +1,193 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// same as line_minimization_03 but use three points cubic fit. + +#include + +#include + +#include +#include + +#include "../tests.h" +using namespace dealii; + +/* + * MWE in Maxima + + +Case 1 (compare to Table I) Function 5.1 Figure 3: + +b:2; +func(x):=-x/(x^2+b); +gfunc(x):=''(diff(func(x),x)); +mu:0.001; +eta : 0.1; +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 ); +pline(x):=func(0)+x*mu*gfunc(0); +plot2d([func(x),pline(x),w(x)], [x,0,16]); + +Case 2, Function 5.2 Figure 4 (compared to Table II we have more iterations, up +to 8 vs 16): + +b:0.004; +func(x):=(x+b)^5-2*(x+b)^4; +gfunc(x):=''(diff(func(x),x)); +mu:0.1; +eta : 0.1; +w1(x) := func(x) - func(0) - x * mu * gfunc(0); +w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0)); +w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 ); +pline(x):=func(0)+x*mu*gfunc(0); +plot2d([func(x),pline(x),w(x)], [x,0,2]); + +Case 3-5 Function 5.4 Figure 6 (Table IV-V-VI): + +b1:0.001; +b2:0.001; + +b1:0.01; +b2:0.001; + +b1:0.001; +b2:0.01; + +g(x):=(1+x^2)^(1/2)-x; +func(x):=g(b1)*((1-x)^2+b2^2)^(1/2) + g(b2)*(x^2+b1^2)^(1/2); +gfunc(x):=''(diff(func(x),x)); +plot2d([func(x)], [x,0,1]); + + */ + + +void +test() +{ + const std::vector values = {{1e-3, 1e-1, 1e+1, 1e+3}}; + double f0, g0, fi, gi; + + { + deallog << "Table 1:" << std::endl; + const double b = 2; + auto func = [&](const double x) { + const double f = -x / (x * x + b); + const double g = 2. * x * x / std::pow(x * x + 2., 2) - 1. / (x * x + 2.); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit_three_points, + a1, + 0.1, + 0.001); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + + { + deallog << "Table 2:" << std::endl; + const double b = 0.004; + auto func = [&](const double x) { + const double f = std::pow(x + b, 5) - 2. * std::pow(x + b, 4); + const double g = 5. * std::pow(x + b, 4) - 8. * std::pow(x + b, 3); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit_three_points, + a1, + 0.100001, + 0.1); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + + { + const std::vector> params = { + {{0.001, 0.001}, {0.01, 0.001}, {0.001, 0.01}}}; + + unsigned int ind = 4; + for (auto p : params) + { + deallog << "Table " << ind++ << ":" << std::endl; + const double b1 = p.first; + const double b2 = p.second; + + const double gb1 = std::sqrt(1. + b1 * b1) - b1; + const double gb2 = std::sqrt(1. + b2 * b2) - b2; + + auto func = [&](const double x) { + const double f = gb1 * std::sqrt(std::pow(1. - x, 2) + b2 * b2) + + gb2 * std::sqrt(x * x + b1 * b1); + const double g = + gb2 * x / sqrt(x * x + b1 * b1) - + gb1 * (1. - x) / std::sqrt(std::pow(1 - x, 2) + b2 * b2); + return std::make_pair(f, g); + }; + + const auto fg0 = func(0); + + for (auto a1 : values) + { + const auto res = LineMinimization::line_search( + func, + fg0.first, + fg0.second, + LineMinimization::poly_fit_three_points, + a1, + 0.00100001, + 0.001); + + const auto fgi = func(res.first); + deallog << res.second << " " << res.first << " " << fgi.second + << std::endl; + } + } + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/line_minimization_03b.output b/tests/optimization/line_minimization_03b.output new file mode 100644 index 0000000000..43e82b9536 --- /dev/null +++ b/tests/optimization/line_minimization_03b.output @@ -0,0 +1,25 @@ +DEAL::Table 1: +DEAL::8 1.49141 0.0125702 +DEAL::6 1.49995 0.0138338 +DEAL::1 10.0000 0.00941945 +DEAL::5 25.6769 0.00150302 +DEAL::Table 2: +DEAL::15 1.59600 2.94690e-10 +DEAL::19 1.59600 9.73245e-10 +DEAL::11 1.59600 -2.83316e-08 +DEAL::19 1.59600 -7.43938e-12 +DEAL::Table 4: +DEAL::3 0.0910000 -5.97089e-05 +DEAL::1 0.100000 -4.93296e-05 +DEAL::3 0.431189 -1.14275e-06 +DEAL::5 0.693314 4.27147e-06 +DEAL::Table 5: +DEAL::6 0.0741039 -2.26077e-05 +DEAL::3 0.0754771 0.000296756 +DEAL::11 0.0734116 -0.000190322 +DEAL::14 0.0723147 -0.000465817 +DEAL::Table 6: +DEAL::13 0.924531 -0.000294927 +DEAL::11 0.925274 -0.000124134 +DEAL::6 0.926593 0.000191539 +DEAL::10 0.928406 0.000653661 diff --git a/tests/optimization/quadratic_fit.cc b/tests/optimization/quadratic_fit.cc new file mode 100644 index 0000000000..8183bb8c78 --- /dev/null +++ b/tests/optimization/quadratic_fit.cc @@ -0,0 +1,64 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//--------------------------------------------------------------- + +// check minimization of the quadratic fit based on f(x1), f(x2) and f'(x1) + +#include + +#include + +#include +#include + +#include "../tests.h" +using namespace dealii; + + +void +test() +{ + // test 1: + { + auto f = [](double x) { + return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x; + }; + auto g = [](double x) { + return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1; + }; + + const double x1 = 10; + const double x2 = 17; + const double f1 = f(x1); + const double f2 = f(x2); + const double g1 = g(x1); + const double g2 = g(x2); + const double res = *LineMinimization::quadratic_fit(x1, f1, g1, x2, f2); + deallog << x1 << " " << x2 << std::endl + << f1 << " " << f2 << std::endl + << g1 << " " << g2 << std::endl + << res << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile, /*do not print job id*/ false); + deallog.depth_console(0); + + test(); +} diff --git a/tests/optimization/quadratic_fit.output b/tests/optimization/quadratic_fit.output new file mode 100644 index 0000000000..8537015bd5 --- /dev/null +++ b/tests/optimization/quadratic_fit.output @@ -0,0 +1,4 @@ +DEAL::10.0000 17.0000 +DEAL::-9999.00 -14737.3 +DEAL::-1999.90 2312.10 +DEAL::15.2907