14#include <dune/common/float_cmp.hh>
15#include <dune/functions/common/signature.hh>
17#include <spdlog/spdlog.h>
19#include <Eigen/Sparse>
46 double maxtime = std::numeric_limits<double>::infinity();
55 double Delta_bar = std::numeric_limits<double>::infinity();
68 static_assert(std::copy_constructible<UF>,
69 " The update function should be copy constructable. If it is a lambda wrap it in a std::function");
78 template <
typename UF2>
85 template <
typename IDBC2>
95 typename UF = utils::UpdateDefault,
typename IDBCF = utils::IDBCForceDefault>
106template <
typename F,
typename TRConfig>
107requires traits::isSpecializationNonTypeAndTypes<TrustRegionConfig, std::remove_cvref_t<TRConfig>>::value
109 static constexpr PreConditioner preConditioner = std::remove_cvref_t<TRConfig>::preConditionerType;
110 using UF = std::remove_cvref_t<TRConfig>::UpdateFunction;
111 using IDBCF = std::remove_cvref_t<TRConfig>::IDBCForceFunction;
112 static_assert(std::remove_cvref_t<F>::nDerivatives == 2,
113 "The number of derivatives in the DifferentiableFunction have to be exactly 2.");
114 auto solver = std::make_shared<TrustRegion<F, preConditioner, UF, IDBCF>>(
115 f, std::forward<TRConfig>(config).updateFunction, std::forward<TRConfig>(config).idbcForceFunction);
117 solver->setup(config.parameters);
183template <
typename F, PreConditioner preConditioner,
typename UF,
typename IDBCF>
206 template <
typename UF2 = UF,
typename IDBCF2 = IDBCF>
208 : energyFunction_{f},
211 static_assert(std::same_as<IDBCForceFunction, utils::IDBCForceDefault>,
212 "Trust Region Method is not implemented to handle inhomogeneous Dirichlet BCs.");
220 settings_ = settings;
221 assert(settings_.
rho_prime < 0.25 &&
"options.rho_prime must be strictly smaller than 1/4.");
222 assert(settings_.
Delta_bar > 0 &&
"options.Delta_bar must be positive.");
224 "options.Delta0 must be positive and smaller than Delta_bar.");
243 auto& energyF = energyFunction_;
244 auto gradientF = derivative(energyF);
245 auto hessianF = derivative(gradientF);
246 decltype(
auto) e = energyF(x);
247 decltype(
auto) g = gradientF(x);
248 decltype(
auto) h = hessianF(x);
252 truncatedConjugateGradient_.analyzePattern(h);
255 truncatedConjugateGradient_.analyzePattern(h);
259 " | iter | inner_i | rho | energy | energy_p | energy_inc | norm(g) | Delta | norm(corr) | "
261 spdlog::info(
"{:-^143}",
"-");
262 while (not stoppingCriterion(e)) {
267 while (eta_.dot(h * eta_) > innerInfo_.
Delta * innerInfo_.
Delta)
274 solveInnerProblem(g, h);
285 const Eigen::VectorXd Hg = h * g;
286 const auto g_Hg = g.dot(Hg);
290 tauC = std::min(Dune::power(stats_.
gradNorm, 3) / (innerInfo_.
Delta * g_Hg), 1.0);
293 const Eigen::VectorXd etaC = -tauC * innerInfo_.
Delta / stats_.
gradNorm * g;
294 const Eigen::VectorXd HetaC = -tauC * innerInfo_.
Delta / stats_.
gradNorm * Hg;
296 const double mdle = stats_.
energy + g.dot(eta_) + .5 * Heta_.dot(eta_);
297 const double mdlec = stats_.
energy + g.dot(etaC) + .5 * HetaC.dot(etaC);
298 if (mdlec < mdle && stats_.
outerIter == 0) {
310 updateFunction_(x, eta_);
319 auto rhoden = -eta_.dot(g + 0.5 * Heta_);
324 const auto rhoReg = std::max(1.0, abs(stats_.
energy)) * eps_ * settings_.
rho_reg;
325 rhonum = rhonum + rhoReg;
326 rhoden = rhoden + rhoReg;
328 const bool modelDecreased = rhoden > 0.0;
333 stats_.
rho = rhonum / rhoden;
334 stats_.
rho = stats_.
rho < 0.0 ? -1.0 : stats_.
rho;
342 if (stats_.
rho < 1e-4 || not modelDecreased || std::isnan(stats_.
rho) || not energyDecreased) {
344 innerInfo_.
Delta /= 4.0;
349 spdlog::info(
" +++ Detected many consecutive TR- (radius decreases).");
350 spdlog::info(
" +++ Consider decreasing options.Delta_bar by an order of magnitude.");
361 spdlog::info(
" +++ Detected many consecutive TR+ (radius increases)");
362 spdlog::info(
" +++ Consider increasing options.Delta_bar by an order of magnitude");
370 if (modelDecreased && stats_.
rho > settings_.
rho_prime && energyDecreased) {
373 "Energy function increased by {} (step size: {}). Since this is small we accept the step and hope for "
374 "convergence of the gradient norm.",
385 innerInfo_.
Delta /= 2;
408 updateFunction_(x, -eta_);
427 return solverInformation;
434 auto&
energy() {
return energyFunction_; }
440 const auto&
energy()
const {
return energyFunction_; }
446 auto residual()
const {
return derivative(energyFunction_); }
462 constexpr auto make_optional_reference(T& value) {
463 return std::make_optional<std::reference_wrapper<const T>>(std::cref(value));
467 requires(not std::is_lvalue_reference_v<T>)
468 constexpr T make_optional_reference(T&& value) {
472 void logState()
const {
474 "{:>3s} {:>3s} {:>6d} {:>9d} {:>6.2f} {:>9.2e} {:>9.2e} {:>11.2e} {:>9.2e} {:>9.2e} {:>11.2e} "
481 void logFinalState() {
482 spdlog::info(
"{:>3s} {:>3s} {:>6d} {:>9d} {: ^6} {: ^9} {: ^9} {: ^11} {:>9.2e} {: ^9} {: ^11} {:<73}",
487 bool stoppingCriterion(
const auto&
energy) {
488 std::ostringstream stream;
492 spdlog::info(
"CONVERGENCE: Energy: {:1.16e} norm(gradient): {:1.16e}",
energy, stats_.
gradNorm);
493 stream <<
"Gradient norm tolerance reached; options.tolerance = " << settings_.
grad_tol;
501 spdlog::info(
"CONVERGENCE: Energy: {:1.16e} norm(correction): {:1.16e}",
energy, stats_.
etaNorm);
502 stream <<
"Displacement norm tolerance reached; = " << settings_.
corr_tol <<
"." << std::endl;
512 stream <<
"Max time exceeded; options.maxtime = " << settings_.
maxtime <<
".";
521 stream <<
"Max iteration count reached; options.maxiter = " << settings_.
maxIter <<
".";
529 void solveInnerProblem(
const auto& g,
const auto& h) {
530 truncatedConjugateGradient_.
setInfo(innerInfo_);
532 truncatedConjugateGradient_.factorize(h);
537 while (truncatedConjugateGradient_.info() != Eigen::Success) {
538 choleskyInitialShift_ *= 5;
539 truncatedConjugateGradient_.preconditioner().setInitialShift(choleskyInitialShift_);
540 truncatedConjugateGradient_.factorize(h);
542 DUNE_THROW(Dune::MathError,
"Factorization of preconditioner failed!");
545 if (truncatedConjugateGradient_.info() == Eigen::Success)
546 choleskyInitialShift_ = 1e-3;
548 eta_ = truncatedConjugateGradient_.solveWithGuess(-g, eta_);
549 innerInfo_ = truncatedConjugateGradient_.
getInfo();
556 std::remove_cvref_t<CorrectionType> eta_;
557 std::remove_cvref_t<CorrectionType> Heta_;
560 double choleskyInitialShift_ = 1e-3;
563 static constexpr double eps_ = 0.0001220703125;
564 std::array<std::string, 6> tcg_stop_reason_{
565 {
"negative curvature",
"exceeded trust region",
"reached target residual-kappa (linear)",
566 "reached target residual-theta (superlinear)",
"maximum inner iterations",
"model increased"}
569 using PreConditionerType =
572 typename Eigen::DiagonalPreconditioner<std::decay_t<EnergyType>>,
573 typename Eigen::IncompleteCholesky<std::decay_t<EnergyType>>>>;
575 truncatedConjugateGradient_;
589 typename UF = utils::UpdateDefault,
typename IDBCF = utils::IDBCForceDefault>
590auto makeTrustRegion(
const F& f, UF&& updateFunction = {}, IDBCF&& idbcForceFunction = {}) {
591 return std::make_shared<TrustRegion<F, preConditioner, UF, IDBCF>>(f, updateFunction, idbcForceFunction);
595 typename UF2 = utils::UpdateDefault,
typename IDBCF2 = utils::IDBCForceDefault>
596TrustRegion(
const F& f, UF2&& updateFunction = {}, IDBCF2&& idbcForceFunction = {})
597 ->
TrustRegion<F, preConditioner, std::remove_cvref_t<UF2>, std::remove_cvref_t<IDBCF2>>;
Helper for the autodiff library.
Contains stl-like type traits.
Collection of fallback default functions.
Implementation of the observer design pattern with broadcasters.
Enums for observer messages.
Implementation of the solver information returned by the nonlinear solvers.
Base for all nonlinear solvers.
Definition of TruncatedConjugateGradient class for solving linear systems using truncated conjugate g...
auto norm(const Eigen::MatrixBase< Derived > &v)
Adding free norm function to Eigen types.
Definition: linearalgebrahelper.hh:258
Definition: assemblermanipulatorbuildingblocks.hh:22
auto makeTrustRegion(const F &f, UF &&updateFunction={}, IDBCF &&idbcForceFunction={})
Creates an instance of the TrustRegion solver.
Definition: trustregion.hh:590
::value auto createNonlinearSolver(NRConfig &&config, F &&f)
Function to create a NewtonRaphson solver instance.
Definition: newtonraphson.hh:106
NonLinearSolverMessages
Enum class defining non-linear solver-related messages.
Definition: broadcastermessages.hh:22
TrustRegion(const F &f, UF2 &&updateFunction={}, IDBCF2 &&idbcForceFunction={}) -> TrustRegion< F, preConditioner, std::remove_cvref_t< UF2 >, std::remove_cvref_t< IDBCF2 > >
StopReason
Enumeration of reasons for stopping the TrustRegion solver.
Definition: trustregion.hh:126
@ maximumIterationsReached
@ correctionNormTolReached
PreConditioner
Enumeration of available preconditioners for the trust region solver.
Definition: trustregion.hh:37
Scalar Delta
Definition: truncatedconjugategradient.hh:38
TCGStopReason stop_tCG
Definition: truncatedconjugategradient.hh:37
Eigen::Index numInnerIter
Definition: truncatedconjugategradient.hh:43
Iterative solver for solving linear systems using the truncated conjugate gradient method.
Definition: truncatedconjugategradient.hh:196
TCGInfo< typename MatrixType::RealScalar > getInfo()
Get information about the truncated conjugate gradient algorithm.
Definition: truncatedconjugategradient.hh:227
void setInfo(TCGInfo< typename MatrixType::RealScalar > alginfo)
Set information about the truncated conjugate gradient algorithm.
Definition: truncatedconjugategradient.hh:233
Base for all nonlinear solvers. Defines the message interface that can be broadcasted to listeners.
Definition: nonlinearsolverbase.hh:24
State for nonlinear solvers.
Definition: nonlinearsolverstate.hh:24
const Domain & domain
Definition: nonlinearsolverstate.hh:28
Information about the result of a non-linear solver.
Definition: solverinfos.hh:21
double correctionNorm
Definition: solverinfos.hh:30
int iterations
Definition: solverinfos.hh:31
double residualNorm
Definition: solverinfos.hh:29
bool success
Definition: solverinfos.hh:28
Definition: trustregion.hh:44
double grad_tol
Gradient tolerance.
Definition: trustregion.hh:50
double Delta_bar
Maximum trust region radius.
Definition: trustregion.hh:55
double rho_reg
Regularization value for rho.
Definition: trustregion.hh:54
double Delta0
Initial trust region radius.
Definition: trustregion.hh:56
bool useRand
Flag for using random correction predictor.
Definition: trustregion.hh:53
int verbosity
Verbosity level.
Definition: trustregion.hh:45
int debug
Debugging flag.
Definition: trustregion.hh:49
int minIter
Minimum number of iterations.
Definition: trustregion.hh:47
double maxtime
Maximum allowable time for solving.
Definition: trustregion.hh:46
double corr_tol
Correction tolerance.
Definition: trustregion.hh:51
int maxIter
Maximum number of iterations.
Definition: trustregion.hh:48
double rho_prime
Rho prime value.
Definition: trustregion.hh:52
Definition: trustregion.hh:67
UF UpdateFunction
Definition: trustregion.hh:71
UF updateFunction
Definition: trustregion.hh:76
auto rebindUpdateFunction(UF2 &&updateFunction) const
Definition: trustregion.hh:79
TRSettings parameters
Definition: trustregion.hh:74
auto rebindIDBCForceFunction(IDBC2 &&idbcForceFunction) const
Definition: trustregion.hh:86
IDBCF idbcForceFunction
Definition: trustregion.hh:77
static constexpr PreConditioner preConditionerType
Definition: trustregion.hh:75
IDBCF IDBCForceFunction
Definition: trustregion.hh:72
Trust Region solver for non-linear optimization problems.
Definition: trustregion.hh:185
TrustRegion(const F &f, UF2 &&updateFunction={}, IDBCF2 &&idbcForceFunction={})
Constructs a TrustRegion solver instance.
Definition: trustregion.hh:207
auto & energy()
Access the energy function.
Definition: trustregion.hh:434
typename FTraits::template Range< 0 > EnergyType
Type of the scalar cost.
Definition: trustregion.hh:196
const IDBCForceFunction & idbcForceFunction() const
Access the force function calculating internal forces due to inhomogeneous Dirichlet BCs.
Definition: trustregion.hh:458
typename FTraits::Domain Domain
Type of the parameter vector of.
Definition: trustregion.hh:190
auto residual() const
Access the residual.
Definition: trustregion.hh:446
IDBCF IDBCForceFunction
Type representing the force function to handle inhomogeneous Dirichlet BCs.
Definition: trustregion.hh:194
HessianType JacobianType
Definition: trustregion.hh:199
typename FTraits::template Range< 1 > GradientType
Type of the gradient vector.
Definition: trustregion.hh:197
F DifferentiableFunction
Type of function to minimize.
Definition: trustregion.hh:193
typename FTraits::template Range< 1 > CorrectionType
Type of the correction of x += deltaX.
Definition: trustregion.hh:191
void setup(const Settings &settings)
Sets up the TrustRegion solver with the provided settings and checks feasibility.
Definition: trustregion.hh:219
TRSettings Settings
Type of the settings for the TrustRegion solver.
Definition: trustregion.hh:187
const UpdateFunction & updateFunction() const
Access the update function.
Definition: trustregion.hh:452
typename F::Traits FTraits
Definition: trustregion.hh:189
NonLinearSolverInformation solve(Domain &x, double stepSize=0.0)
Solves the nonlinear optimization problem using the TrustRegion algorithm.
Definition: trustregion.hh:233
UF UpdateFunction
Type of the update function.
Definition: trustregion.hh:192
const auto & energy() const
Access the energy function.
Definition: trustregion.hh:440
typename FTraits::template Range< 2 > HessianType
Type of the Hessian matrix.
Definition: trustregion.hh:198
Additional information about the TrustRegion algorithm.
Definition: trustregion.hh:139
std::string cauchystr
Used Cauchy step string.
Definition: trustregion.hh:147
StopReason stop
Stopping reason.
Definition: trustregion.hh:150
std::string randomPredictionString
Random prediction string.
Definition: trustregion.hh:146
bool acceptProposal
Flag indicating whether the proposal is accepted.
Definition: trustregion.hh:148
std::string reasonString
String describing the stopping reason.
Definition: trustregion.hh:151
int consecutive_TRminus
Consecutive trust region decreases.
Definition: trustregion.hh:141
bool used_cauchy
Flag indicating whether Cauchy point was used.
Definition: trustregion.hh:149
int consecutive_Rejected
Consecutive rejected proposals.
Definition: trustregion.hh:142
std::string trstr
Trust region change string (TR+, TR-).
Definition: trustregion.hh:144
std::string stopReasonString
String describing the stopping reason.
Definition: trustregion.hh:143
std::string accstr
Acceptance string (REJ, acc).
Definition: trustregion.hh:145
int consecutive_TRplus
Consecutive trust region increases.
Definition: trustregion.hh:140
Information about the TrustRegion solver.
Definition: trustregion.hh:159
double rho
Definition: trustregion.hh:165
double etaNorm
Definition: trustregion.hh:161
double energyProposal
Definition: trustregion.hh:164
double energy
Definition: trustregion.hh:163
double gradNorm
Definition: trustregion.hh:160
int outerIter
Definition: trustregion.hh:166
double time
Definition: trustregion.hh:162
int innerIterSum
Definition: trustregion.hh:167
void notify(NonLinearSolverMessages message, const State &data)
This calls all the registered functions.
Definition: broadcaster.hh:61
Default struct used to represent that no inhomogeneous Dirichlet BCs are present.
Definition: defaultfunctions.hh:49
Default functor for updating operations.
Definition: defaultfunctions.hh:65