13#include <dune/common/float_cmp.hh>
15#include <spdlog/spdlog.h>
17#include <Eigen/Sparse>
46 double maxtime = std::numeric_limits<double>::infinity();
55 double Delta_bar = std::numeric_limits<double>::infinity();
121 typename UF = utils::UpdateDefault>
125 using ValueType =
typename NLO::template ParameterValue<0>;
132 using ScalarType = std::remove_cvref_t<typename NLO::template FunctionReturnType<0>>;
135 using MatrixType = std::remove_cvref_t<typename NLO::template FunctionReturnType<2>>;
144 updateFunction_{updateFunction},
146 eta_.setZero(gradient().size());
147 Heta_.setZero(gradient().size());
148 truncatedConjugateGradient_.analyzePattern(hessian());
157 assert(options_.
rho_prime < 0.25 &&
"options.rho_prime must be strictly smaller than 1/4.");
158 assert(options_.
Delta_bar > 0 &&
"options.Delta_bar must be positive.");
160 "options.Delta0 must be positive and smaller than Delta_bar.");
174 template <
typename SolutionType = NoPredictor>
175 requires std::is_same_v<SolutionType, NoPredictor> || std::is_convertible_v<SolutionType, CorrectionType>
181 NonLinearSolverInformation solverInformation;
187 if constexpr (not std::is_same_v<SolutionType, NoPredictor>)
188 updateFunction(x, dxPredictor);
189 truncatedConjugateGradient_.analyzePattern(hessian());
193 " | iter | inner_i | rho | energy | energy_p | energy_inc | norm(g) | Delta | norm(corr) | "
195 spdlog::info(
"{:-^143}",
"-");
196 while (not stoppingCriterion()) {
201 while (eta_.dot(hessian() * eta_) > innerInfo_.
Delta * innerInfo_.
Delta)
212 Heta_ = hessian() * eta_;
219 const Eigen::VectorXd Hg = hessian() * gradient();
220 const auto g_Hg = (gradient().dot(Hg));
224 tauC = std::min(Dune::power(stats_.
gradNorm, 3) / (innerInfo_.
Delta * g_Hg), 1.0);
227 const Eigen::VectorXd etaC = -tauC * innerInfo_.
Delta / stats_.
gradNorm * gradient();
228 const Eigen::VectorXd HetaC = -tauC * innerInfo_.
Delta / stats_.
gradNorm * Hg;
230 const double mdle = stats_.
energy + gradient().dot(eta_) + .5 * Heta_.dot(eta_);
231 const double mdlec = stats_.
energy + gradient().dot(etaC) + .5 * HetaC.dot(etaC);
232 if (mdlec < mdle && stats_.
outerIter == 0) {
244 updateFunction_(x, eta_);
253 auto rhoden = -eta_.dot(gradient() + 0.5 * Heta_);
258 const auto rhoReg = std::max(1.0, abs(stats_.
energy)) * eps_ * options_.
rho_reg;
259 rhonum = rhonum + rhoReg;
260 rhoden = rhoden + rhoReg;
262 const bool modelDecreased = rhoden > 0.0;
267 stats_.
rho = rhonum / rhoden;
268 stats_.
rho = stats_.
rho < 0.0 ? -1.0 : stats_.
rho;
276 if (stats_.
rho < 1e-4 || not modelDecreased || std::isnan(stats_.
rho) || not energyDecreased) {
278 innerInfo_.
Delta /= 4.0;
283 spdlog::info(
" +++ Detected many consecutive TR- (radius decreases).");
284 spdlog::info(
" +++ Consider decreasing options.Delta_bar by an order of magnitude.");
295 spdlog::info(
" +++ Detected many consecutive TR+ (radius increases)");
296 spdlog::info(
" +++ Consider increasing options.Delta_bar by an order of magnitude");
304 if (modelDecreased && stats_.
rho > options_.
rho_prime && energyDecreased) {
307 "Energy function increased by {} (step size: {}). Since this is small we accept the step and hope for "
308 "convergence of the gradient norm.",
319 innerInfo_.
Delta /= 2;
334 nonLinearOperator_.updateAll();
343 nonLinearOperator_.updateAll();
344 stats_.
gradNorm = gradient().norm();
350 solverInformation.success =
353 solverInformation.iterations = stats_.
outerIter;
354 solverInformation.residualNorm = stats_.
gradNorm;
355 if (solverInformation.success)
357 return solverInformation;
366 void logState()
const {
368 "{:>3s} {:>3s} {:>6d} {:>9d} {:>6.2f} {:>9.2e} {:>9.2e} {:>11.2e} {:>9.2e} {:>9.2e} {:>11.2e} "
375 void logFinalState() {
376 spdlog::info(
"{:>3s} {:>3s} {:>6d} {:>9d} {: ^6} {: ^9} {: ^9} {: ^11} {:>9.2e} {: ^9} {: ^11} {:<73}",
385 bool stoppingCriterion() {
386 std::ostringstream stream;
390 spdlog::info(
"CONVERGENCE: Energy: {:1.16e} norm(gradient): {:1.16e}",
nonLinearOperator().value(),
392 stream <<
"Gradient norm tolerance reached; options.tolerance = " << options_.
grad_tol;
400 spdlog::info(
"CONVERGENCE: Energy: {:1.16e} norm(correction): {:1.16e}",
nonLinearOperator().value(),
402 stream <<
"Displacement norm tolerance reached; = " << options_.
corr_tol <<
"." << std::endl;
412 stream <<
"Max time exceeded; options.maxtime = " << options_.
maxtime <<
".";
421 stream <<
"Max iteration count reached; options.maxiter = " << options_.
maxIter <<
".";
429 void solveInnerProblem() {
430 truncatedConjugateGradient_.
setInfo(innerInfo_);
432 truncatedConjugateGradient_.factorize(hessian());
437 while (truncatedConjugateGradient_.info() != Eigen::Success) {
438 choleskyInitialShift_ *= 5;
439 truncatedConjugateGradient_.preconditioner().setInitialShift(choleskyInitialShift_);
440 truncatedConjugateGradient_.factorize(hessian());
442 DUNE_THROW(Dune::MathError,
"Factorization of preconditioner failed!");
445 if (truncatedConjugateGradient_.info() == Eigen::Success)
446 choleskyInitialShift_ = 1e-3;
448 eta_ = truncatedConjugateGradient_.solveWithGuess(-gradient(), eta_);
449 innerInfo_ = truncatedConjugateGradient_.
getInfo();
452 NLO nonLinearOperator_;
454 typename NLO::template ParameterValue<0> xOld_;
457 TrustRegionSettings options_;
459 double choleskyInitialShift_ = 1e-3;
462 static constexpr double eps_ = 0.0001220703125;
463 std::array<std::string, 6> tcg_stop_reason_{
464 {
"negative curvature",
"exceeded trust region",
"reached target residual-kappa (linear)",
465 "reached target residual-theta (superlinear)",
"maximum inner iterations",
"model increased"}
468 using PreConditionerType =
471 typename Eigen::DiagonalPreconditioner<ScalarType>,
472 typename Eigen::IncompleteCholesky<ScalarType>>>;
474 truncatedConjugateGradient_;
488 typename UF = utils::UpdateDefault>
490 return std::make_shared<TrustRegion<NLO, preConditioner, UF>>(nonLinearOperator, updateFunction);
Contains stl-like type traits.
Collection of fallback default functions.
Helper for the autodiff library.
Implementation of the Newton-Raphson method for solving nonlinear equations.
Definition of TruncatedConjugateGradient class for solving linear systems using truncated conjugate g...
Enums for observer messages.
Implementation of the observer design pattern.
auto norm(const Eigen::MatrixBase< Derived > &v)
Adding free norm function to Eigen types.
Definition: linearalgebrahelper.hh:259
Definition: simpleassemblers.hh:22
auto makeTrustRegion(const NLO &nonLinearOperator, UF &&updateFunction={})
Creates an instance of the TrustRegion solver.
Definition: trustregion.hh:489
StopReason
Enumeration of reasons for stopping the TrustRegion solver.
Definition: trustregion.hh:64
@ maximumIterationsReached
@ correctionNormTolReached
PreConditioner
Enumeration of available preconditioners for the trust region solver.
Definition: trustregion.hh:34
Scalar Delta
Definition: truncatedconjugategradient.hh:38
TCGStopReason stop_tCG
Definition: truncatedconjugategradient.hh:37
Eigen::Index numInnerIter
Definition: truncatedconjugategradient.hh:43
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
Information about the result of a non-linear solver.
Definition: solverinfos.hh:19
Configuration settings for the TrustRegion solver.
Definition: trustregion.hh:44
double rho_reg
Regularization value for rho.
Definition: trustregion.hh:54
int minIter
Minimum number of iterations.
Definition: trustregion.hh:47
int verbosity
Verbosity level.
Definition: trustregion.hh:45
int maxIter
Maximum number of iterations.
Definition: trustregion.hh:48
double grad_tol
Gradient tolerance.
Definition: trustregion.hh:50
double Delta0
Initial trust region radius.
Definition: trustregion.hh:56
int debug
Debugging flag.
Definition: trustregion.hh:49
double Delta_bar
Maximum trust region radius.
Definition: trustregion.hh:55
bool useRand
Flag for using random correction predictor.
Definition: trustregion.hh:53
double rho_prime
Rho prime value.
Definition: trustregion.hh:52
double maxtime
Maximum allowable time for solving.
Definition: trustregion.hh:46
double corr_tol
Correction tolerance.
Definition: trustregion.hh:51
Additional information about the TrustRegion algorithm.
Definition: trustregion.hh:77
std::string cauchystr
Used Cauchy step string.
Definition: trustregion.hh:85
StopReason stop
Stopping reason.
Definition: trustregion.hh:88
std::string randomPredictionString
Random prediction string.
Definition: trustregion.hh:84
bool acceptProposal
Flag indicating whether the proposal is accepted.
Definition: trustregion.hh:86
std::string reasonString
String describing the stopping reason.
Definition: trustregion.hh:89
int consecutive_TRminus
Consecutive trust region decreases.
Definition: trustregion.hh:79
bool used_cauchy
Flag indicating whether Cauchy point was used.
Definition: trustregion.hh:87
int consecutive_Rejected
Consecutive rejected proposals.
Definition: trustregion.hh:80
std::string trstr
Trust region change string (TR+, TR-).
Definition: trustregion.hh:82
std::string stopReasonString
String describing the stopping reason.
Definition: trustregion.hh:81
std::string accstr
Acceptance string (REJ, acc).
Definition: trustregion.hh:83
int consecutive_TRplus
Consecutive trust region increases.
Definition: trustregion.hh:78
Information about the TrustRegion solver.
Definition: trustregion.hh:97
double rho
Definition: trustregion.hh:103
double etaNorm
Definition: trustregion.hh:99
double energyProposal
Definition: trustregion.hh:102
double energy
Definition: trustregion.hh:101
double gradNorm
Definition: trustregion.hh:98
int outerIter
Definition: trustregion.hh:104
double time
Definition: trustregion.hh:100
int innerIterSum
Definition: trustregion.hh:105
Trust Region solver for non-linear optimization problems.
Definition: trustregion.hh:123
typename NLO::DerivativeType CorrectionType
Type of the correction of x += deltaX.
Definition: trustregion.hh:127
std::remove_cvref_t< typename NLO::template FunctionReturnType< 0 > > ScalarType
Definition: trustregion.hh:133
typename NLO::template ParameterValue< 0 > ValueType
Definition: trustregion.hh:126
void setup(const TrustRegionSettings &settings)
Sets up the TrustRegion solver with the provided settings and checks feasibility.
Definition: trustregion.hh:155
auto & nonLinearOperator()
Access the nonlinear operator.
Definition: trustregion.hh:363
std::remove_cvref_t< typename NLO::template FunctionReturnType< 2 > > MatrixType
Type of the Hessian.
Definition: trustregion.hh:135
NLO NonLinearOperator
Type of the non-linear operator.
Definition: trustregion.hh:130
UF UpdateFunction
Type of the update function.
Definition: trustregion.hh:128
TrustRegion(const NLO &nonLinearOperator, UF updateFunction={})
Constructs a TrustRegion solver instance.
Definition: trustregion.hh:142
NonLinearSolverInformation solve(const SolutionType &dxPredictor=NoPredictor{})
Solves the nonlinear optimization problem using the TrustRegion algorithm.
Definition: trustregion.hh:176
Generic observable interface for the Observer design pattern. See for a description of the design pa...
Definition: observer.hh:129
void notify(NonLinearSolverMessages message)
Notify observers about a specific message type.