version 0.4.4
trustregion.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
9#pragma once
10
11#include <iosfwd>
12#include <type_traits>
13
14#include <dune/common/float_cmp.hh>
15#include <dune/functions/common/signature.hh>
16
17#include <spdlog/spdlog.h>
18
19#include <Eigen/Sparse>
20
29
30namespace Ikarus {
31
37{
41};
42
44{
45 int verbosity = 5;
46 double maxtime = std::numeric_limits<double>::infinity();
47 int minIter = 3;
48 int maxIter = 1000;
49 int debug = 0;
50 double grad_tol = 1e-6;
51 double corr_tol = 1e-6;
52 double rho_prime = 0.01;
53 bool useRand = false;
54 double rho_reg = 1e6;
55 double Delta_bar = std::numeric_limits<double>::infinity();
56 double Delta0 = 10;
57};
58
64template <PreConditioner preConditioner = PreConditioner::IncompleteCholesky, typename UF = utils::UpdateDefault,
65 typename IDBCF = utils::IDBCForceDefault>
67{
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");
70
71 using UpdateFunction = UF;
72 using IDBCForceFunction = IDBCF;
73
75 static constexpr PreConditioner preConditionerType = preConditioner;
78 template <typename UF2>
81 .updateFunction = std::forward<UF2>(updateFunction),
82 .idbcForceFunction = idbcForceFunction};
83 return settings;
84 }
85 template <typename IDBC2>
88 .updateFunction = updateFunction,
89 .idbcForceFunction = std::forward<IDBC2>(idbcForceFunction)};
90 return settings;
91 }
92};
93
94template <typename F, PreConditioner preConditioner = PreConditioner::IncompleteCholesky,
95 typename UF = utils::UpdateDefault, typename IDBCF = utils::IDBCForceDefault>
96class TrustRegion;
97
106template <typename F, typename TRConfig>
107requires traits::isSpecializationNonTypeAndTypes<TrustRegionConfig, std::remove_cvref_t<TRConfig>>::value
108auto createNonlinearSolver(TRConfig&& config, F&& f) {
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);
116
117 solver->setup(config.parameters);
118 return solver;
119}
120
125enum class StopReason
126{
132};
133
139{
143 std::string stopReasonString;
144 std::string trstr;
145 std::string accstr;
147 std::string cauchystr = " ";
149 bool used_cauchy = false;
151 std::string reasonString;
152};
153
158struct Stats
159{
160 double gradNorm{1};
161 double etaNorm{1};
162 double time{0};
163 double energy{0};
164 double energyProposal{0};
165 double rho{0};
166 int outerIter{0};
168};
169
183template <typename F, PreConditioner preConditioner, typename UF, typename IDBCF>
185{
186public:
188
189 using FTraits = typename F::Traits;
190 using Domain = typename FTraits::Domain;
191 using CorrectionType = typename FTraits::template Range<1>;
192 using UpdateFunction = UF;
194 using IDBCForceFunction = IDBCF;
195
196 using EnergyType = typename FTraits::template Range<0>;
197 using GradientType = typename FTraits::template Range<1>;
198 using HessianType = typename FTraits::template Range<2>;
206 template <typename UF2 = UF, typename IDBCF2 = IDBCF>
207 explicit TrustRegion(const F& f, UF2&& updateFunction = {}, IDBCF2&& idbcForceFunction = {})
208 : energyFunction_{f},
209 updateFunction_{std::forward<UF2>(updateFunction)},
210 idbcForceFunction_{std::forward<IDBCF2>(idbcForceFunction)} {
211 static_assert(std::same_as<IDBCForceFunction, utils::IDBCForceDefault>,
212 "Trust Region Method is not implemented to handle inhomogeneous Dirichlet BCs.");
213 }
214
219 void setup(const Settings& settings) {
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.");
223 assert(settings_.Delta0 > 0 && settings_.Delta0 < settings_.Delta_bar &&
224 "options.Delta0 must be positive and smaller than Delta_bar.");
225 }
226
233 [[nodiscard]] NonLinearSolverInformation solve(Domain& x, double stepSize = 0.0) {
234 using enum NonLinearSolverMessages;
235
236 NonLinearSolverInformation solverInformation;
237 auto state = typename TrustRegion::State{.domain = x, .correction = eta_, .information = solverInformation};
238
239 this->notify(INIT, state);
240 stats_ = Stats{};
241 info_ = AlgoInfo{};
242
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);
249
250 eta_.resizeLike(g);
251 Heta_.resizeLike(g);
252 truncatedConjugateGradient_.analyzePattern(h);
253 stats_.energy = e;
254 stats_.gradNorm = norm(g);
255 truncatedConjugateGradient_.analyzePattern(h);
256
257 innerInfo_.Delta = settings_.Delta0;
258 spdlog::info(
259 " | iter | inner_i | rho | energy | energy_p | energy_inc | norm(g) | Delta | norm(corr) | "
260 "InnerBreakReason");
261 spdlog::info("{:-^143}", "-");
262 while (not stoppingCriterion(e)) {
263 this->notify(ITERATION_STARTED, state);
264 if (settings_.useRand) {
265 if (stats_.outerIter == 0) {
266 eta_.setRandom();
267 while (eta_.dot(h * eta_) > innerInfo_.Delta * innerInfo_.Delta)
268 eta_ *= eps_; // eps is sqrt(sqrt(maschine-precision))
269 } else
270 eta_.setZero();
271 } else
272 eta_.setZero();
273
274 solveInnerProblem(g, h);
275 stats_.innerIterSum += innerInfo_.numInnerIter;
276
277 info_.stopReasonString = tcg_stop_reason_[static_cast<int>(innerInfo_.stop_tCG)];
278 Heta_ = h * eta_;
279 if (settings_.useRand and stats_.outerIter == 0) {
280 info_.used_cauchy = false;
281 info_.randomPredictionString = " Used Random correction predictor";
282 info_.cauchystr = " ";
283 double tauC;
284 // Check the curvature
285 const Eigen::VectorXd Hg = h * g;
286 const auto g_Hg = g.dot(Hg);
287 if (g_Hg <= 0)
288 tauC = 1;
289 else
290 tauC = std::min(Dune::power(stats_.gradNorm, 3) / (innerInfo_.Delta * g_Hg), 1.0);
291
292 // generate the Cauchy point.
293 const Eigen::VectorXd etaC = -tauC * innerInfo_.Delta / stats_.gradNorm * g;
294 const Eigen::VectorXd HetaC = -tauC * innerInfo_.Delta / stats_.gradNorm * Hg;
295
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) {
299 eta_ = etaC;
300 Heta_ = HetaC;
301 info_.used_cauchy = true;
302
303 info_.cauchystr = " USED CAUCHY POINT!";
304 }
305 } else
306 info_.cauchystr = " ";
307
308 stats_.etaNorm = eta_.norm();
309
310 updateFunction_(x, eta_);
311
312 // Calculate energy of our proposed update step
313 e = energyF(x);
314 stats_.energyProposal = e;
315
316 // Will we accept the proposal or not?
317 // Check the performance of the quadratic model against the actual energy.
318 auto rhonum = stats_.energy - stats_.energyProposal;
319 auto rhoden = -eta_.dot(g + 0.5 * Heta_);
320
321 /* Close to convergence the proposed energy and the real energy almost coincide.
322 * Therefore, the performance check of our model becomes ill-conditioned
323 * The regularisation fixes this */
324 const auto rhoReg = std::max(1.0, abs(stats_.energy)) * eps_ * settings_.rho_reg;
325 rhonum = rhonum + rhoReg;
326 rhoden = rhoden + rhoReg;
327
328 const bool modelDecreased = rhoden > 0.0;
329
330 if (!modelDecreased)
331 info_.stopReasonString.append(", model did not decrease");
332
333 stats_.rho = rhonum / rhoden;
334 stats_.rho = stats_.rho < 0.0 ? -1.0 : stats_.rho; // move rho to the domain [-1.0,inf]
335
336 info_.trstr = " ";
337
338 // measure if energy decreased
339 const bool energyDecreased = Dune::FloatCmp::ge(stats_.energy - stats_.energyProposal, -1e-12);
340
341 // If the model behaves badly or if the energy increased we reduce the trust region size
342 if (stats_.rho < 1e-4 || not modelDecreased || std::isnan(stats_.rho) || not energyDecreased) {
343 info_.trstr = "TR-";
344 innerInfo_.Delta /= 4.0;
345 info_.consecutive_TRplus = 0;
346 info_.consecutive_TRminus++;
347 if (info_.consecutive_TRminus >= 5 && settings_.verbosity >= 1) {
348 info_.consecutive_TRminus = -std::numeric_limits<int>::infinity();
349 spdlog::info(" +++ Detected many consecutive TR- (radius decreases).");
350 spdlog::info(" +++ Consider decreasing options.Delta_bar by an order of magnitude.");
351 }
352
353 } else if (stats_.rho > 0.99 && (innerInfo_.stop_tCG == Eigen::TCGStopReason::negativeCurvature ||
355 info_.trstr = "TR+";
356 innerInfo_.Delta = std::min(3.5 * innerInfo_.Delta, settings_.Delta_bar);
357 info_.consecutive_TRminus = 0;
358 info_.consecutive_TRplus++;
359 if (info_.consecutive_TRplus >= 5 && settings_.verbosity >= 1) {
360 info_.consecutive_TRplus = -std::numeric_limits<int>::infinity();
361 spdlog::info(" +++ Detected many consecutive TR+ (radius increases)");
362 spdlog::info(" +++ Consider increasing options.Delta_bar by an order of magnitude");
363
364 } else {
365 info_.consecutive_TRplus = 0;
366 info_.consecutive_TRminus = 0;
367 }
368 }
369
370 if (modelDecreased && stats_.rho > settings_.rho_prime && energyDecreased) {
371 if (stats_.energyProposal > stats_.energy)
372 spdlog::info(
373 "Energy function increased by {} (step size: {}). Since this is small we accept the step and hope for "
374 "convergence of the gradient norm.",
375 stats_.energyProposal - stats_.energy, stats_.etaNorm);
376
377 info_.acceptProposal = true;
378 info_.accstr = "acc";
379 info_.consecutive_Rejected = 0;
380 } else {
381 info_.acceptProposal = false;
382 info_.accstr = "REJ";
383
384 if (info_.consecutive_Rejected >= 5)
385 innerInfo_.Delta /= 2;
386 else
387 innerInfo_.Delta = std::min(innerInfo_.Delta, stats_.etaNorm / 2.0);
388 ++info_.consecutive_Rejected;
389 }
390
391 stats_.outerIter++;
392
393 if (settings_.verbosity == 1)
394 logState();
395
396 info_.randomPredictionString = "";
397
398 solverInformation.correctionNorm = stats_.etaNorm;
399 solverInformation.residualNorm = stats_.gradNorm;
400 this->notify(CORRECTION_UPDATED, state);
401
402 if (info_.acceptProposal) {
403 stats_.energy = stats_.energyProposal;
404 this->notify(CORRECTIONNORM_UPDATED, state);
405 this->notify(RESIDUALNORM_UPDATED, state);
406 this->notify(SOLUTION_CHANGED, state);
407 } else {
408 updateFunction_(x, -eta_);
409 eta_.setZero();
410 }
411 e = energyF(x);
412 g = gradientF(x);
413 h = hessianF(x);
414 stats_.gradNorm = g.norm();
416 }
417 spdlog::info("{}", info_.reasonString);
418 spdlog::info("Total iterations: {} Total CG Iterations: {}", stats_.outerIter, stats_.innerIterSum);
419
420 solverInformation.success =
422
423 solverInformation.iterations = stats_.outerIter;
424 solverInformation.residualNorm = stats_.gradNorm;
425 if (solverInformation.success)
427 return solverInformation;
428 }
429
434 auto& energy() { return energyFunction_; }
435
440 const auto& energy() const { return energyFunction_; }
441
446 auto residual() const { return derivative(energyFunction_); }
447
452 const UpdateFunction& updateFunction() const { return updateFunction_; }
453
458 const IDBCForceFunction& idbcForceFunction() const { return idbcForceFunction_; }
459
460private:
461 template <class T>
462 constexpr auto make_optional_reference(T& value) {
463 return std::make_optional<std::reference_wrapper<const T>>(std::cref(value));
464 }
465
466 template <class T>
467 requires(not std::is_lvalue_reference_v<T>)
468 constexpr T make_optional_reference(T&& value) {
469 return value;
470 }
471
472 void logState() const {
473 spdlog::info(
474 "{:>3s} {:>3s} {:>6d} {:>9d} {:>6.2f} {:>9.2e} {:>9.2e} {:>11.2e} {:>9.2e} {:>9.2e} {:>11.2e} "
475 "{:<73}",
476 info_.accstr, info_.trstr, stats_.outerIter, innerInfo_.numInnerIter, stats_.rho, stats_.energy,
477 stats_.energyProposal, stats_.energyProposal - stats_.energy, stats_.gradNorm, innerInfo_.Delta, stats_.etaNorm,
479 }
480
481 void logFinalState() {
482 spdlog::info("{:>3s} {:>3s} {:>6d} {:>9d} {: ^6} {: ^9} {: ^9} {: ^11} {:>9.2e} {: ^9} {: ^11} {:<73}",
483 info_.accstr, info_.trstr, stats_.outerIter, innerInfo_.numInnerIter, " ", " ", " ", " ",
484 stats_.gradNorm, " ", " ", info_.stopReasonString + info_.cauchystr + info_.randomPredictionString);
485 }
486
487 bool stoppingCriterion(const auto& energy) {
488 std::ostringstream stream;
490 if (stats_.gradNorm < settings_.grad_tol && stats_.outerIter != 0) {
491 logFinalState();
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;
494
495 info_.reasonString = stream.str();
496
498 return true;
499 } else if (stats_.etaNorm < settings_.corr_tol && stats_.outerIter != 0) {
500 logFinalState();
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;
503
504 info_.reasonString = stream.str();
506 return true;
507 }
508
510 if (stats_.time >= settings_.maxtime) {
511 logFinalState();
512 stream << "Max time exceeded; options.maxtime = " << settings_.maxtime << ".";
513 info_.reasonString = stream.str();
515 return true;
516 }
517
519 if (stats_.outerIter >= settings_.maxIter) {
520 logFinalState();
521 stream << "Max iteration count reached; options.maxiter = " << settings_.maxIter << ".";
522 info_.reasonString = stream.str();
524 return true;
525 }
526 return false;
527 }
528
529 void solveInnerProblem(const auto& g, const auto& h) {
530 truncatedConjugateGradient_.setInfo(innerInfo_);
531 int attempts = 0;
532 truncatedConjugateGradient_.factorize(h);
533 // If the preconditioner is IncompleteCholesky the factorization may fail if we have negative diagonal entries and
534 // the initial shift is too small. Therefore, if the factorization fails we increase the initial shift by a factor
535 // of 5.
536 if constexpr (preConditioner == PreConditioner::IncompleteCholesky) {
537 while (truncatedConjugateGradient_.info() != Eigen::Success) {
538 choleskyInitialShift_ *= 5;
539 truncatedConjugateGradient_.preconditioner().setInitialShift(choleskyInitialShift_);
540 truncatedConjugateGradient_.factorize(h);
541 if (attempts > 5)
542 DUNE_THROW(Dune::MathError, "Factorization of preconditioner failed!");
543 ++attempts;
544 }
545 if (truncatedConjugateGradient_.info() == Eigen::Success)
546 choleskyInitialShift_ = 1e-3;
547 }
548 eta_ = truncatedConjugateGradient_.solveWithGuess(-g, eta_);
549 innerInfo_ = truncatedConjugateGradient_.getInfo();
550 }
551
552 F energyFunction_;
553
554 UpdateFunction updateFunction_;
555 IDBCForceFunction idbcForceFunction_;
556 std::remove_cvref_t<CorrectionType> eta_;
557 std::remove_cvref_t<CorrectionType> Heta_;
558 Settings settings_;
559 AlgoInfo info_;
560 double choleskyInitialShift_ = 1e-3;
561 Eigen::TCGInfo<double> innerInfo_;
562 Stats stats_;
563 static constexpr double eps_ = 0.0001220703125; // 0.0001220703125 is sqrt(sqrt(maschine-precision))
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"}
567 };
568
569 using PreConditionerType =
570 std::conditional_t<preConditioner == PreConditioner::IdentityPreconditioner, Eigen::IdentityPreconditioner,
571 std::conditional_t<preConditioner == PreConditioner::DiagonalPreconditioner,
572 typename Eigen::DiagonalPreconditioner<std::decay_t<EnergyType>>,
573 typename Eigen::IncompleteCholesky<std::decay_t<EnergyType>>>>;
574 Eigen::TruncatedConjugateGradient<std::decay_t<HessianType>, Eigen::Lower | Eigen::Upper, PreConditionerType>
575 truncatedConjugateGradient_;
576};
577
588template <typename F, PreConditioner preConditioner = PreConditioner::IncompleteCholesky,
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);
592}
593
594template <typename F, PreConditioner preConditioner = PreConditioner::IncompleteCholesky,
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>>;
598
599} // namespace Ikarus
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
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