From: Jean-Paul Pelteret Date: Tue, 18 Sep 2018 09:11:16 +0000 (+0200) Subject: Add code examples for line search. X-Git-Tag: v9.1.0-rc1~696^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5e52252ac520d763fecae9a6f3a6cda36727a97c;p=dealii.git Add code examples for line search. --- diff --git a/include/deal.II/optimization/line_minimization.h b/include/deal.II/optimization/line_minimization.h index f89f51a3d0..494c6be87a 100644 --- a/include/deal.II/optimization/line_minimization.h +++ b/include/deal.II/optimization/line_minimization.h @@ -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( + * func, + * fg0.first, fg0.second, + * LineMinimization::poly_fit, + * 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 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 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 solution_total_trial = + * get_solution_total(newton_update_trial); + * + * // Recompute the linear system based on the trial newton update + * Vector system_rhs (...); + * SparseMatrix tangent_matrix (...); + * assemble_linear_system( + * tangent_matrix, system_rhs, solution_total_trial); + * Vector 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 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( + * ls_minimization_function, + * res_0.first, res_0.second, + * LineMinimization::poly_fit, + * 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.