Reference documentation for deal.II version 9.4.1
\(\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
Functions
LineMinimization Namespace Reference

Functions

template<typename NumberType >
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)
 
template<typename NumberType >
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)
 
template<typename NumberType >
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)
 
template<typename NumberType >
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)
 
template<typename NumberType >
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)
 
template<typename NumberType >
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)
 

Detailed Description

A namespace for various algorithms related to minimization a over line.

Function Documentation

◆ quadratic_fit()

template<typename NumberType >
std_cxx17::optional< NumberType > LineMinimization::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 the gradient \(g(x\_low)\), return the local minimizer of the quadratic interpolation function.

The return type is optional to fit with similar functions that may not have a solution for given parameters.

◆ cubic_fit()

template<typename NumberType >
std_cxx17::optional< NumberType > LineMinimization::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 
)

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 value).

The return type is optional as the real-valued solution might not exist.

◆ cubic_fit_three_points()

template<typename NumberType >
std_cxx17::optional< NumberType > LineMinimization::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 
)

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.

◆ poly_fit()

template<typename NumberType >
NumberType LineMinimization::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 
)

Return the minimizer of a polynomial using function values f_low , f_hi , and f_rec[0] at three points x_low , x_hi , and x_rec[0] as well as the derivatives at two points g_low and g_hi. The returned point should be within the bounds bounds .

This function will first try to perform a cubic_fit(). If its unsuccessful, or if the minimum is not within the provided bounds, a quadratic_fit() will be performed. The function will fallback to a bisection method if quadratic_fit() fails as well.

◆ poly_fit_three_points()

template<typename NumberType >
NumberType LineMinimization::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 
)

Same as poly_fit(), but performing a cubic fit with three points (see cubic_fit_three_points() ).

◆ line_search()

template<typename NumberType >
std::pair< NumberType, unsigned int > LineMinimization::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 
)

Perform a line search in \((0,max]\) with strong Wolfe conditions

\[ f(\alpha) \le f(0) + \alpha \mu f'(0) \\ |f'(\alpha)| \le \eta |f'(0)| \]

using the one dimensional function func in conjunction with a function interpolate to choose a new point from the interval based on the function values and derivatives at its ends. The parameter a1 is a trial estimate of the first step. Interpolation can be done using either poly_fit() or poly_fit_three_points(), or any other function that has a similar signature.

The function implements Algorithms 2.6.2 and 2.6.4 on pages 34-35 in

@book{Fletcher2013,
title = {Practical methods of optimization},
publisher = {John Wiley \& Sons},
year = {2013},
author = {Fletcher, Roger},
isbn = {978-0-471-49463-8},
doi = {10.1002/9781118723203},
}

These are minor variations of Algorithms 3.5 and 3.6 on pages 60-61 in

@book{Nocedal2006,
title = {Numerical Optimization},
publisher = {Springer New York},
year = {2006},
author = {Jorge Nocedal and S. Wright},
address = {233 Spring Street, New York, NY 10013, USA},
doi = {10.1007/978-0-387-40065-5},
}

It consists of a bracketing phase and a zoom phase, where interpolate is used.

Two examples of use might be as follows: In the first example, we wish to find the minimum of the function \(100 * x^4 + (1-x)^2\). To find the approximate solution using line search with a polynomial fit to the curve one would perform the following steps:

auto func = [](const double x)
{
const double f = 100. * std::pow(x, 4) + std::pow(1. - x, 2); // Value
const double g = 400. * std::pow(x, 3) - 2. * (1. - x); // Gradient
return std::make_pair(f, g);
};
const auto fg0 = func(0);
const auto res = LineMinimization::line_search<double>(
func,
fg0.first, fg0.second,
LineMinimization::poly_fit<double>,
0.1, 0.1, 0.01, 100, 20);
const double approx_solution = res.first;
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition: vectorization.h:5526

In the second example, we wish to perform line search in the context of a non-linear finite element problem. What follows below is a non-optimized implementation of the back-tracking algorithm, which may be useful when the load-step size is too large. The following illustrates the basic steps necessary to utilize the scheme within the context of a global nonlinear solver:

// Solve some incremental linear system
const Vector<double> newton_update = solver_linear_system(...);
// Now we check to see if the suggested Newton update is a good one.
// First we define what it means to perform linesearch in the context of
// this incremental nonlinear finite element problem.
auto ls_minimization_function = [&](const double step_size)
{
// Scale the full Newton update by the proposed line search step size.
Vector<double> newton_update_trial(newton_update);
newton_update_trial *= step_size;
// Ensure that the Dirichlet constraints are correctly applied,
// irrespective of the step size
constraints.distribute(newton_update_trial);
// Now add the constribution from the previously accepted solution
// history.
const Vector<double> solution_total_trial =
get_solution_total(newton_update_trial);
// Recompute the linear system based on the trial newton update
Vector<double> system_rhs (...);
SparseMatrix<double> tangent_matrix (...);
assemble_linear_system(
tangent_matrix, system_rhs, solution_total_trial);
Vector<double> residual_trial (system_rhs);
residual_trial *= -1.0; // Residual = -RHS
// Negelect the constrained entries in the consideration
// of the function (value and gradient) to be minimized.
constraints.set_zero(residual_trial);
// Here we compute the function value according to the text given in
// section 5.1.4 of Wriggers, P., "Nonlinear finite element methods",
// 2008.
// The function value correspeonds to equ. 5.11 on p159.
const double f = 0.5 * (residual_trial * residual_trial); // Value
// However, the corresponding gradient given in eq 5.14 is wrong. The
// suggested result
// const double g = -(residual_0*residual_trial);
// should actually be
// g = G(V + alpha*delta)*[ K(V + alpha*delta)*delta.
Vector<double> tmp;
tmp.reinit(newton_update);
tangent_matrix.vmult(tmp, newton_update);
const double g = tmp * residual_trial; // Gradient
return std::make_pair(f, g);
};
// Next we can write a function to determine if taking the full Newton
// step is a good idea or not (i.e. if it offers good convergence
// characterisics). This function calls the one we defined above,
// and actually only performs the line search if an early exit
// criterion is not met.
auto perform_linesearch = [&]()
{
const auto res_0 = ls_minimization_function(0.0);
Assert(res_0.second < 0.0,
ExcMessage("Gradient should be negative. Current value: " +
std::to_string(res_0.second)));
const auto res_1 = ls_minimization_function(1.0);
// Check to see if the minimum lies in the interval [0,1] through the
// values of the gradients at the limit points.
// If it does not, then the full step is accepted. This is discussed by
// Wriggers in the paragraph after equ. 5.14.
if (res_0.second * res_1.second > 0.0)
return 1.0;
// The values for eta, mu are chosen such that more strict convergence
// conditions are enforced.
// They should be adjusted according to the problem requirements.
const double a1 = 1.0;
const double eta = 0.5;
const double mu = 0.49;
const double a_max = 1.25;
const double max_evals = 20;
const auto res = LineMinimization::line_search<double>(
ls_minimization_function,
res_0.first, res_0.second,
LineMinimization::poly_fit<double>,
a1, eta, mu, a_max, max_evals));
return res.first; // Final stepsize
};
// Finally, we can perform the line search and adjust the Newton update
// accordingly.
const double linesearch_step_size = perform_linesearch();
if (linesearch_step_size != 1.0)
{
newton_update *= linesearch_step_size;
constraints.distribute(newton_update);
}
Vector
Definition: vector.h:109
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
Vector::reinit
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Parameters
funcA one dimensional function which returns value and derivative at the given point.
f0The function value at the origin.
g0The function derivative at the origin.
interpolateA 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.
a1Initial trial step for the bracketing phase.
etaA parameter in the second Wolfe condition (curvature condition).
muA parameter in the first Wolfe condition (sufficient decrease).
a_maxThe maximum allowed step size.
max_evaluationsThe maximum allowed number of function evaluations.
debug_outputA flag to output extra debug information into the deallog static object.
Returns
The function returns the step size and the number of times function func was called.
Generated by doxygen 1.9.5