]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add code examples for line search. 7182/head
authorJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 18 Sep 2018 09:11:16 +0000 (11:11 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 18 Sep 2018 09:27:58 +0000 (11:27 +0200)
include/deal.II/optimization/line_minimization.h

index f89f51a3d0f188ff92b01057eb7aeb8b4943efb4..494c6be87ac68fa5fb8177611993dae89b1c7608 100644 (file)
@@ -173,6 +173,135 @@ namespace LineMinimization
    * @endcode
    * It consists of a bracketing phase and a zoom phase, where @p 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:
+   *
+   * @code
+   *   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;
+   * @endcode
+   *
+   * 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:
+   *
+   * @code
+   *   // 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
+   *     // <code>const double g = -(residual_0*residual_trial);</code>
+   *     // 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);
+   *   }
+   * @endcode
+   *
+   *
    * @param func A one dimensional function which returns value and derivative
    * at the given point.
    * @param f0 The function value at the origin.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.