ROL
ROL_TypeU_TrustRegionAlgorithm_Def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Rapid Optimization Library (ROL) Package
4//
5// Copyright 2014 NTESS and the ROL contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ROL_TRUSTREGIONALGORITHM_U_DEF_H
11#define ROL_TRUSTREGIONALGORITHM_U_DEF_H
12
14
15namespace ROL {
16namespace TypeU {
17
18template<typename Real>
20 const Ptr<Secant<Real>> &secant )
21 : Algorithm<Real>(), esec_(SECANT_USERDEFINED) {
22 // Set status test
23 status_->reset();
24 status_->add(makePtr<StatusTest<Real>>(parlist));
25
26 // Trust-Region Parameters
27 ParameterList &slist = parlist.sublist("Step");
28 ParameterList &trlist = slist.sublist("Trust Region");
29 state_->searchSize = trlist.get("Initial Radius", static_cast<Real>(-1));
30 delMax_ = trlist.get("Maximum Radius", ROL_INF<Real>());
31 eta0_ = trlist.get("Step Acceptance Threshold", static_cast<Real>(0.05));
32 eta1_ = trlist.get("Radius Shrinking Threshold", static_cast<Real>(0.05));
33 eta2_ = trlist.get("Radius Growing Threshold", static_cast<Real>(0.9));
34 gamma0_ = trlist.get("Radius Shrinking Rate (Negative rho)", static_cast<Real>(0.0625));
35 gamma1_ = trlist.get("Radius Shrinking Rate (Positive rho)", static_cast<Real>(0.25));
36 gamma2_ = trlist.get("Radius Growing Rate", static_cast<Real>(2.5));
37 TRsafe_ = trlist.get("Safeguard Size", static_cast<Real>(100.0));
38 eps_ = TRsafe_*ROL_EPSILON<Real>();
39 // Nonmonotone Information
40 NMstorage_ = trlist.get("Nonmonotone Storage Limit", 0);
41 useNM_ = (NMstorage_ <= 0 ? false : true);
42 // Inexactness Information
43 ParameterList &glist = parlist.sublist("General");
44 useInexact_.clear();
45 useInexact_.push_back(glist.get("Inexact Objective Function", false));
46 useInexact_.push_back(glist.get("Inexact Gradient", false));
47 useInexact_.push_back(glist.get("Inexact Hessian-Times-A-Vector", false));
48 // Trust-Region Inexactness Parameters
49 ParameterList &ilist = trlist.sublist("Inexact").sublist("Gradient");
50 scale0_ = ilist.get("Tolerance Scaling", static_cast<Real>(0.1));
51 scale1_ = ilist.get("Relative Tolerance", static_cast<Real>(2));
52 // Inexact Function Evaluation Information
53 ParameterList &vlist = trlist.sublist("Inexact").sublist("Value");
54 scale_ = vlist.get("Tolerance Scaling", static_cast<Real>(1.e-1));
55 omega_ = vlist.get("Exponent", static_cast<Real>(0.9));
56 force_ = vlist.get("Forcing Sequence Initial Value", static_cast<Real>(1.0));
57 updateIter_ = vlist.get("Forcing Sequence Update Frequency", static_cast<int>(10));
58 forceFactor_ = vlist.get("Forcing Sequence Reduction Factor", static_cast<Real>(0.1));
59 // Initialize Trust Region Subproblem Solver Object
60 std::string solverName = trlist.get("Subproblem Solver", "Dogleg");
61 etr_ = StringToETrustRegionU(solverName);
62 solver_ = TrustRegionUFactory<Real>(parlist);
63 verbosity_ = glist.get("Output Level", 0);
64 // Secant Information
65 useSecantPrecond_ = glist.sublist("Secant").get("Use as Preconditioner", false);
66 useSecantHessVec_ = glist.sublist("Secant").get("Use as Hessian", false);
67 if (secant == nullPtr) {
68 std::string secantType = glist.sublist("Secant").get("Type","Limited-Memory BFGS");
69 esec_ = StringToESecant(secantType);
70 }
71 // Initialize trust region model
72 model_ = makePtr<TrustRegionModel_U<Real>>(parlist,secant);
73 printHeader_ = verbosity_ > 2;
74}
75
76template<typename Real>
77void TrustRegionAlgorithm<Real>::initialize( const Vector<Real> &x,
78 const Vector<Real> &g,
79 Vector<Real> &Bg,
80 Objective<Real> &obj,
81 std::ostream &outStream) {
82 const Real zero(0);
83 // Initialize data
84 Algorithm<Real>::initialize(x,g);
85 solver_->initialize(x,g);
86 model_->initialize(x,g);
87 // Update approximate gradient and approximate objective function.
88 Real ftol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>();
89 obj.update(x,UpdateType::Initial,state_->iter);
90 state_->value = obj.value(x,ftol);
91 state_->nfval++;
92 state_->snorm = ROL_INF<Real>();
93 state_->gnorm = ROL_INF<Real>();
94 Real Delta = state_->searchSize;
95 if (Delta <= zero) state_->searchSize = 1e2*x.norm();
96 computeGradient(x,obj,true);
97 // Check if inverse Hessian is implemented for dogleg methods
98 model_->validate(obj,x,g,etr_);
99 // Compute initial trust region radius if desired.
100 if ( Delta <= zero ) {
101 int nfval = 0;
102 state_->searchSize
103 = TRUtils::initialRadius<Real>(nfval,x,*state_->gradientVec,Bg,
104 state_->value,state_->gnorm,gtol_,obj,*model_,delMax_,
105 outStream,(verbosity_>1));
106 state_->nfval += nfval;
107 }
108}
109
110template<typename Real>
111Real TrustRegionAlgorithm<Real>::computeValue( const Vector<Real> &x,
112 Objective<Real> &obj,
113 Real pRed) {
114 const Real one(1);
115 Real tol(std::sqrt(ROL_EPSILON<Real>())), fval(0);
116 if ( useInexact_[0] ) {
117 if ( !(state_->iter%updateIter_) && (state_->iter != 0) ) {
118 force_ *= forceFactor_;
119 }
120 Real eta = static_cast<Real>(0.999)*std::min(eta1_,one-eta2_);
121 tol = scale_*std::pow(eta*std::min(pRed,force_),one/omega_);
122 state_->value = obj.value(*state_->iterateVec,tol);
123 state_->nfval++;
124 }
125 // Evaluate objective function at new iterate
126 obj.update(x,UpdateType::Trial);
127 fval = obj.value(x,tol);
128 state_->nfval++;
129 return fval;
130}
131
132template<typename Real>
133void TrustRegionAlgorithm<Real>::computeGradient(const Vector<Real> &x,
134 Objective<Real> &obj,
135 bool accept) {
136 if ( useInexact_[1] ) {
137 Real gtol0 = scale0_*state_->searchSize;
138 if (accept) gtol_ = gtol0 + static_cast<Real>(1);
139 else gtol0 = scale0_*std::min(state_->gnorm,state_->searchSize);
140 while ( gtol_ > gtol0 ) {
141 gtol_ = gtol0;
142 obj.gradient(*state_->gradientVec,x,gtol_); state_->ngrad++;
143 state_->gnorm = state_->gradientVec->norm();
144 gtol0 = scale0_*std::min(state_->gnorm,state_->searchSize);
145 }
146 }
147 else {
148 if (accept) {
149 gtol_ = std::sqrt(ROL_EPSILON<Real>());
150 obj.gradient(*state_->gradientVec,x,gtol_); state_->ngrad++;
151 state_->gnorm = state_->gradientVec->norm();
152 }
153 }
154}
155
156template<typename Real>
157void TrustRegionAlgorithm<Real>::run( Vector<Real> &x,
158 const Vector<Real> &g,
159 Objective<Real> &obj,
160 std::ostream &outStream ) {
161 const Real zero(0);
162 // Initialize trust-region data
163 Real ftrial(0), pRed(0), rho(0);
164 Ptr<Vector<Real>> gvec = g.clone();
165 initialize(x,g,*gvec,obj,outStream);
166 // Initialize nonmonotone data
167 Real rhoNM(0), sigmac(0), sigmar(0);
168 Real fr(state_->value), fc(state_->value), fmin(state_->value);
169 TRUtils::ETRFlag TRflagNM;
170 int L(0);
171
172 // Output
173 if (verbosity_ > 0) writeOutput(outStream,true);
174
175 while (status_->check(*state_)) {
176 // Build trust-region model
177 model_->setData(obj,x,*state_->gradientVec,gtol_);
178 // Minimize trust-region model over trust-region constraint
179 pRed = zero;
180 SPflag_ = 0; SPiter_ = 0;
181 solver_->solve(*state_->stepVec,state_->snorm,pRed,SPflag_,SPiter_,
182 state_->searchSize,*model_);
183 // Compute trial objective function value
184 x.plus(*state_->stepVec);
185 ftrial = computeValue(x,obj,pRed);
186 // Compute ratio of actual and predicted reduction
187 TRflag_ = TRUtils::SUCCESS;
188 TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
189 if (useNM_) {
190 TRUtils::analyzeRatio<Real>(rhoNM,TRflagNM,fr,ftrial,pRed+sigmar,eps_,outStream,verbosity_>1);
191 TRflag_ = (rho < rhoNM ? TRflagNM : TRflag_);
192 rho = (rho < rhoNM ? rhoNM : rho );
193 }
194 // Update algorithm state
195 state_->iter++;
196 // Accept/reject step and update trust region radius
197 if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS)
198 || (TRflag_ >= 2)) { // Step Rejected
199 x.set(*state_->iterateVec);
200 obj.update(x,UpdateType::Revert,state_->iter);
201 if (rho < zero && TRflag_ != TRUtils::TRNAN) {
202 // Negative reduction, interpolate to find new trust-region radius
203 state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
204 state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
205 outStream,verbosity_>1);
206 }
207 else { // Shrink trust-region radius
208 state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
209 }
210 computeGradient(x,obj,false);
211 }
212 else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
213 || (TRflag_ == TRUtils::POSPREDNEG)) { // Step Accepted
214 state_->iterateVec->set(x);
215 state_->value = ftrial;
216 obj.update(x,UpdateType::Accept,state_->iter);
217 if (useNM_) {
218 sigmac += pRed; sigmar += pRed;
219 if (ftrial < fmin) { fmin = ftrial; fc = fmin; sigmac = zero; L = 0; }
220 else {
221 L++;
222 if (ftrial > fc) { fc = ftrial; sigmac = zero; }
223 if (L == NMstorage_) { fr = fc; sigmar = sigmac; }
224 }
225 }
226 // Increase trust-region radius
227 if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
228 // Compute gradient at new iterate
229 gvec->set(*state_->gradientVec);
230 computeGradient(x,obj,true);
231 // Update secant information in trust-region model
232 model_->update(x,*state_->stepVec,*gvec,*state_->gradientVec,
233 state_->snorm,state_->iter);
234 }
235 // Update Output
236 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
237 }
238 if (verbosity_ > 0) Algorithm<Real>::writeExitStatus(outStream);
239}
240
241template<typename Real>
242void TrustRegionAlgorithm<Real>::writeHeader( std::ostream& os ) const {
243 std::ios_base::fmtflags osFlags(os.flags());
244 if(verbosity_ > 1) {
245 os << std::string(114,'-') << std::endl;
246 os << "Trust-Region status output definitions" << std::endl << std::endl;
247 os << " iter - Number of iterates (steps taken)" << std::endl;
248 os << " value - Objective function value" << std::endl;
249 os << " gnorm - Norm of the gradient" << std::endl;
250 os << " snorm - Norm of the step (update to optimization vector)" << std::endl;
251 os << " delta - Trust-Region radius" << std::endl;
252 os << " #fval - Number of times the objective function was evaluated" << std::endl;
253 os << " #grad - Number of times the gradient was computed" << std::endl;
254 os << std::endl;
255 os << " tr_flag - Trust-Region flag" << std::endl;
256 for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
257 os << " " << NumberToString(flag) << " - "
258 << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
259 }
260 if( etr_ == TRUSTREGION_U_TRUNCATEDCG ) {
261 os << std::endl;
262 os << " iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
263 os << " flagGC - Trust-Region Truncated CG flag" << std::endl;
264 for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
265 os << " " << NumberToString(flag) << " - "
266 << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
267 }
268 }
269 else if( etr_ == TRUSTREGION_U_SPG ) {
270 os << std::endl;
271 os << " iterCG - Number of spectral projected gradient iterations" << std::endl << std::endl;
272 os << " flagGC - Trust-Region spectral projected gradient flag" << std::endl;
273 }
274 os << std::string(114,'-') << std::endl;
275 }
276 os << " ";
277 os << std::setw(6) << std::left << "iter";
278 os << std::setw(15) << std::left << "value";
279 os << std::setw(15) << std::left << "gnorm";
280 os << std::setw(15) << std::left << "snorm";
281 os << std::setw(15) << std::left << "delta";
282 os << std::setw(10) << std::left << "#fval";
283 os << std::setw(10) << std::left << "#grad";
284 os << std::setw(10) << std::left << "tr_flag";
285 if ( etr_ == TRUSTREGION_U_TRUNCATEDCG ) {
286 os << std::setw(10) << std::left << "iterCG";
287 os << std::setw(10) << std::left << "flagCG";
288 }
289 else if (etr_ == TRUSTREGION_U_SPG) {
290 os << std::setw(10) << std::left << "iterSPG";
291 os << std::setw(10) << std::left << "flagSPG";
292 }
293 os << std::endl;
294 os.flags(osFlags);
295}
296
297template<typename Real>
298void TrustRegionAlgorithm<Real>::writeName( std::ostream& os ) const {
299 std::ios_base::fmtflags osFlags(os.flags());
300 os << std::endl << ETrustRegionUToString(etr_) << " Trust-Region Solver";
301 if ( useSecantPrecond_ || useSecantHessVec_ ) {
302 if ( useSecantPrecond_ && !useSecantHessVec_ ) {
303 os << " with " << ESecantToString(esec_) << " Preconditioning" << std::endl;
304 }
305 else if ( !useSecantPrecond_ && useSecantHessVec_ ) {
306 os << " with " << ESecantToString(esec_) << " Hessian Approximation" << std::endl;
307 }
308 else {
309 os << " with " << ESecantToString(esec_) << " Preconditioning and Hessian Approximation" << std::endl;
310 }
311 }
312 else {
313 os << std::endl;
314 }
315 os.flags(osFlags);
316}
317
318template<typename Real>
319void TrustRegionAlgorithm<Real>::writeOutput(std::ostream& os, bool print_header) const {
320 std::ios_base::fmtflags osFlags(os.flags());
321 os << std::scientific << std::setprecision(6);
322 if ( state_->iter == 0 ) {
323 writeName(os);
324 }
325 if ( print_header ) {
326 writeHeader(os);
327 }
328 if ( state_->iter == 0 ) {
329 os << " ";
330 os << std::setw(6) << std::left << state_->iter;
331 os << std::setw(15) << std::left << state_->value;
332 os << std::setw(15) << std::left << state_->gnorm;
333 os << std::setw(15) << std::left << "---";
334 os << std::setw(15) << std::left << state_->searchSize;
335 os << std::setw(10) << std::left << state_->nfval;
336 os << std::setw(10) << std::left << state_->ngrad;
337 os << std::setw(10) << std::left << "---";
338 if ( etr_ == TRUSTREGION_U_TRUNCATEDCG || etr_ == TRUSTREGION_U_SPG ) {
339 os << std::setw(10) << std::left << "---";
340 os << std::setw(10) << std::left << "---";
341 }
342 os << std::endl;
343 }
344 else {
345 os << " ";
346 os << std::setw(6) << std::left << state_->iter;
347 os << std::setw(15) << std::left << state_->value;
348 os << std::setw(15) << std::left << state_->gnorm;
349 os << std::setw(15) << std::left << state_->snorm;
350 os << std::setw(15) << std::left << state_->searchSize;
351 os << std::setw(10) << std::left << state_->nfval;
352 os << std::setw(10) << std::left << state_->ngrad;
353 os << std::setw(10) << std::left << TRflag_;
354 if ( etr_ == TRUSTREGION_U_TRUNCATEDCG || etr_ == TRUSTREGION_U_SPG ) {
355 os << std::setw(10) << std::left << SPiter_;
356 os << std::setw(10) << std::left << SPflag_;
357 }
358 os << std::endl;
359 }
360 os.flags(osFlags);
361}
362} // namespace TypeU
363} // namespace ROL
364
365#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
virtual void initialize(const Vector< Real > &x)
Initialize temporary variables.
Contains definitions of enums for trust region algorithms.
Provides an interface to run unconstrained optimization algorithms.
TrustRegionAlgorithm(ParameterList &parlist, const Ptr< Secant< Real > > &secant=nullPtr)
std::string NumberToString(T Number)
ETrustRegionU StringToETrustRegionU(std::string s)
ESecant StringToESecant(std::string s)
std::string ETrustRegionUToString(ETrustRegionU tr)
ECGFlag
Enumation of flags used by conjugate gradient methods.
std::string ESecantToString(ESecant tr)
std::string ECGFlagToString(ECGFlag cgf)
@ SECANT_USERDEFINED