LineMinimization Namespace Reference

Functions
template<typename NumberType >
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)
template<typename NumberType >
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)
template<typename NumberType >
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)
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::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::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::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 [87]. These are minor variations of Algorithms 3.5 and 3.6 on pages 60-61 in [165]. 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;

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 contribution 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
// characteristics). 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);
}

Parameters

func	A one dimensional function which returns value and derivative at the given point.
f0	The function value at the origin.
g0	The function derivative at the origin.
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.
a1	Initial trial step for the bracketing phase.
eta	A parameter in the second Wolfe condition (curvature condition).
mu	A parameter in the first Wolfe condition (sufficient decrease).
a_max	The maximum allowed step size.
max_evaluations	The maximum allowed number of function evaluations.
debug_output	A 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.