ROL
ROL_TrustRegion.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_TRUSTREGION_H
45#define ROL_TRUSTREGION_H
46
51#include "ROL_Types.hpp"
56
57namespace ROL {
58
59template<class Real>
61private:
62
63 ROL::Ptr<Vector<Real> > prim_, dual_, xtmp_;
64
66
69 Real pRed_;
71 Real mu0_;
72
73 std::vector<bool> useInexact_;
74
76
79
80 unsigned verbosity_;
81
82 // POST SMOOTHING PARAMETERS
85 Real mu_;
86 Real beta_;
87
88public:
89
90 virtual ~TrustRegion() {}
91
92 // Constructor
93 TrustRegion( ROL::ParameterList &parlist )
94 : pRed_(0), ftol_old_(ROL_OVERFLOW<Real>()), cnt_(0), verbosity_(0) {
95 // Trust-Region Parameters
96 ROL::ParameterList list = parlist.sublist("Step").sublist("Trust Region");
97 TRmodel_ = StringToETrustRegionModel(list.get("Subproblem Model", "Kelley-Sachs"));
98 eta0_ = list.get("Step Acceptance Threshold", static_cast<Real>(0.05));
99 eta1_ = list.get("Radius Shrinking Threshold", static_cast<Real>(0.05));
100 eta2_ = list.get("Radius Growing Threshold", static_cast<Real>(0.9));
101 gamma0_ = list.get("Radius Shrinking Rate (Negative rho)", static_cast<Real>(0.0625));
102 gamma1_ = list.get("Radius Shrinking Rate (Positive rho)", static_cast<Real>(0.25));
103 gamma2_ = list.get("Radius Growing Rate", static_cast<Real>(2.5));
104 mu0_ = list.get("Sufficient Decrease Parameter", static_cast<Real>(1.e-4));
105 TRsafe_ = list.get("Safeguard Size", static_cast<Real>(100.0));
106 eps_ = TRsafe_*ROL_EPSILON<Real>();
107 // General Inexactness Information
108 ROL::ParameterList &glist = parlist.sublist("General");
109 useInexact_.clear();
110 useInexact_.push_back(glist.get("Inexact Objective Function", false));
111 useInexact_.push_back(glist.get("Inexact Gradient", false));
112 useInexact_.push_back(glist.get("Inexact Hessian-Times-A-Vector", false));
113 // Inexact Function Evaluation Information
114 ROL::ParameterList &ilist = list.sublist("Inexact").sublist("Value");
115 scale_ = ilist.get("Tolerance Scaling", static_cast<Real>(1.e-1));
116 omega_ = ilist.get("Exponent", static_cast<Real>(0.9));
117 force_ = ilist.get("Forcing Sequence Initial Value", static_cast<Real>(1.0));
118 updateIter_ = ilist.get("Forcing Sequence Update Frequency", static_cast<int>(10));
119 forceFactor_ = ilist.get("Forcing Sequence Reduction Factor", static_cast<Real>(0.1));
120 // Get verbosity level
121 verbosity_ = glist.get("Print Verbosity", 0);
122 // Post-smoothing parameters
123 max_fval_ = list.sublist("Post-Smoothing").get("Function Evaluation Limit", 20);
124 alpha_init_ = list.sublist("Post-Smoothing").get("Initial Step Size", static_cast<Real>(1));
125 mu_ = list.sublist("Post-Smoothing").get("Tolerance", static_cast<Real>(0.9999));
126 beta_ = list.sublist("Post-Smoothing").get("Rate", static_cast<Real>(0.01));
127 }
128
129 virtual void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
130 prim_ = x.clone();
131 dual_ = g.clone();
132 xtmp_ = x.clone();
133 }
134
135 virtual void update( Vector<Real> &x,
136 Real &fnew,
137 Real &del,
138 int &nfval,
139 int &ngrad,
140 ETrustRegionFlag &flagTR,
141 const Vector<Real> &s,
142 const Real snorm,
143 const Real fold,
144 const Vector<Real> &g,
145 int iter,
146 Objective<Real> &obj,
148 TrustRegionModel<Real> &model ) {
149 Real tol = std::sqrt(ROL_EPSILON<Real>());
150 const Real one(1), zero(0);
151
152 /***************************************************************************************************/
153 // BEGIN INEXACT OBJECTIVE FUNCTION COMPUTATION
154 /***************************************************************************************************/
155 // Update inexact objective function
156 Real fold1 = fold, ftol = tol; // TOL(1.e-2);
157 if ( useInexact_[0] ) {
158 if ( !(cnt_%updateIter_) && (cnt_ != 0) ) {
160 }
161 //const Real oe4(1e4);
162 //Real c = scale_*std::max(TOL,std::min(one,oe4*std::max(pRed_,std::sqrt(ROL_EPSILON<Real>()))));
163 //ftol = c*std::pow(std::min(eta1_,one-eta2_)
164 // *std::min(std::max(pRed_,std::sqrt(ROL_EPSILON<Real>())),force_),one/omega_);
165 //if ( ftol_old_ > ftol || cnt_ == 0 ) {
166 // ftol_old_ = ftol;
167 // fold1 = obj.value(x,ftol_old_);
168 //}
169 //cnt_++;
170 Real eta = static_cast<Real>(0.999)*std::min(eta1_,one-eta2_);
171 ftol = scale_*std::pow(eta*std::min(pRed_,force_),one/omega_);
172 ftol_old_ = ftol;
173 fold1 = obj.value(x,ftol_old_);
174 cnt_++;
175 }
176 // Evaluate objective function at new iterate
177 prim_->set(x); prim_->plus(s);
178 if (bnd.isActivated()) {
179 bnd.project(*prim_);
180 }
181 obj.update(*prim_);
182 fnew = obj.value(*prim_,ftol);
183
184 nfval = 1;
185 Real aRed = fold1 - fnew;
186 /***************************************************************************************************/
187 // FINISH INEXACT OBJECTIVE FUNCTION COMPUTATION
188 /***************************************************************************************************/
189
190 /***************************************************************************************************/
191 // BEGIN COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
192 /***************************************************************************************************/
193 // Modify Actual and Predicted Reduction According to Model
194 model.updateActualReduction(aRed,s);
196
197 if ( verbosity_ > 0 ) {
198 std::cout << std::endl;
199 std::cout << " Computation of actual and predicted reduction" << std::endl;
200 std::cout << " Current objective function value: " << fold1 << std::endl;
201 std::cout << " New objective function value: " << fnew << std::endl;
202 std::cout << " Actual reduction: " << aRed << std::endl;
203 std::cout << " Predicted reduction: " << pRed_ << std::endl;
204 }
205
206 // Compute Ratio of Actual and Predicted Reduction
207 Real EPS = eps_*((one > std::abs(fold1)) ? one : std::abs(fold1));
208 Real aRed_safe = aRed + EPS, pRed_safe = pRed_ + EPS;
209 Real rho(0);
210 if (((std::abs(aRed_safe) < eps_) && (std::abs(pRed_safe) < eps_)) || aRed == pRed_) {
211 rho = one;
213 }
214 else if ( std::isnan(aRed_safe) || std::isnan(pRed_safe) ) {
215 rho = -one;
216 flagTR = TRUSTREGION_FLAG_NAN;
217 }
218 else {
219 rho = aRed_safe/pRed_safe;
220 if (pRed_safe < zero && aRed_safe > zero) {
222 }
223 else if (aRed_safe <= zero && pRed_safe > zero) {
225 }
226 else if (aRed_safe <= zero && pRed_safe < zero) {
228 }
229 else {
231 }
232 }
233
234 if ( verbosity_ > 0 ) {
235 std::cout << " Safeguard: " << eps_ << std::endl;
236 std::cout << " Actual reduction with safeguard: " << aRed_safe << std::endl;
237 std::cout << " Predicted reduction with safeguard: " << pRed_safe << std::endl;
238 std::cout << " Ratio of actual and predicted reduction: " << rho << std::endl;
239 std::cout << " Trust-region flag: " << flagTR << std::endl;
240 }
241 /***************************************************************************************************/
242 // FINISH COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
243 /***************************************************************************************************/
244
245 /***************************************************************************************************/
246 // BEGIN CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
247 /***************************************************************************************************/
248 bool decr = true;
250 if ( rho >= eta0_ && (std::abs(aRed_safe) > eps_) ) {
251 // Compute Criticality Measure || x - P( x - g ) ||
252 prim_->set(x);
253 prim_->axpy(-one,g.dual());
254 bnd.project(*prim_);
255 prim_->scale(-one);
256 prim_->plus(x);
257 Real pgnorm = prim_->norm();
258 // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
259 prim_->set(g.dual());
260 bnd.pruneActive(*prim_,g,x);
261 Real lam = std::min(one, del/prim_->norm());
262 prim_->scale(-lam);
263 prim_->plus(x);
264 bnd.project(*prim_);
265 prim_->scale(-one);
266 prim_->plus(x);
267 pgnorm *= prim_->norm();
268 // Sufficient decrease?
269 decr = ( aRed_safe >= mu0_*pgnorm );
270 flagTR = (!decr ? TRUSTREGION_FLAG_QMINSUFDEC : flagTR);
271
272 if ( verbosity_ > 0 ) {
273 std::cout << " Decrease lower bound (constraints): " << mu0_*pgnorm << std::endl;
274 std::cout << " Trust-region flag (constraints): " << flagTR << std::endl;
275 std::cout << " Is step feasible: " << bnd.isFeasible(x) << std::endl;
276 }
277 }
278 }
279 /***************************************************************************************************/
280 // FINISH CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
281 /***************************************************************************************************/
282
283 /***************************************************************************************************/
284 // BEGIN STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
285 /***************************************************************************************************/
286 if ( verbosity_ > 0 ) {
287 std::cout << " Norm of step: " << snorm << std::endl;
288 std::cout << " Trust-region radius before update: " << del << std::endl;
289 }
290 if ((rho < eta0_ && flagTR == TRUSTREGION_FLAG_SUCCESS) || flagTR >= 2 || !decr ) { // Step Rejected
291 fnew = fold1; // This is a bug if rho < zero...
292 if (rho < zero) { // Negative reduction, interpolate to find new trust-region radius
293 Real gs(0);
294 if ( bnd.isActivated() ) {
295 model.dualTransform(*dual_, *model.getGradient());
296 gs = dual_->dot(s.dual());
297 }
298 else {
299 gs = g.dot(s.dual());
300 }
301 Real modelVal = model.value(s,tol);
302 modelVal += fold1;
303 Real theta = (one-eta2_)*gs/((one-eta2_)*(fold1+gs)+eta2_*modelVal-fnew);
304 del = std::min(gamma1_*std::min(snorm,del),std::max(gamma0_,theta)*del);
305 if ( verbosity_ > 0 ) {
306 std::cout << " Interpolation model value: " << modelVal << std::endl;
307 std::cout << " Interpolation step length: " << theta << std::endl;
308 }
309 }
310 else { // Shrink trust-region radius
311 del = gamma1_*std::min(snorm,del);
312 }
313 obj.update(x,true,iter);
314 }
315 else if ((rho >= eta0_ && flagTR != TRUSTREGION_FLAG_NPOSPREDNEG) ||
316 (flagTR == TRUSTREGION_FLAG_POSPREDNEG)) { // Step Accepted
317 // Perform line search (smoothing) to ensure decrease
319 // Compute new gradient
320 xtmp_->set(x); xtmp_->plus(s);
321 bnd.project(*xtmp_);
322 obj.update(*xtmp_);
323 obj.gradient(*dual_,*xtmp_,tol); // MUST DO SOMETHING HERE WITH TOL
324 ngrad++;
325 // Compute smoothed step
326 Real alpha(1);
327 prim_->set(*xtmp_);
328 prim_->axpy(-alpha/alpha_init_,dual_->dual());
329 bnd.project(*prim_);
330 // Compute new objective value
331 obj.update(*prim_);
332 Real ftmp = obj.value(*prim_,tol); // MUST DO SOMETHING HERE WITH TOL
333 nfval++;
334 // Perform smoothing
335 int cnt = 0;
336 alpha = alpha_init_;
337 while ( (ftmp-fnew) >= mu_*aRed ) {
338 prim_->set(*xtmp_);
339 prim_->axpy(-alpha/alpha_init_,dual_->dual());
340 bnd.project(*prim_);
341 obj.update(*prim_);
342 ftmp = obj.value(*prim_,tol); // MUST DO SOMETHING HERE WITH TOL
343 nfval++;
344 if ( cnt >= max_fval_ ) {
345 break;
346 }
347 alpha *= beta_;
348 cnt++;
349 }
350 // Store objective function and iteration information
351 if (std::isnan(ftmp)) {
352 flagTR = TRUSTREGION_FLAG_NAN;
353 del = gamma1_*std::min(snorm,del);
354 rho = static_cast<Real>(-1);
355 //x.axpy(static_cast<Real>(-1),s);
356 //obj.update(x,true,iter);
357 fnew = fold1;
358 }
359 else {
360 fnew = ftmp;
361 x.set(*prim_);
362 }
363 }
364 else {
365 x.plus(s);
366 }
367 if (rho >= eta2_) { // Increase trust-region radius
368 del = gamma2_*del;
369 }
370 obj.update(x,true,iter);
371 }
372
373 if ( verbosity_ > 0 ) {
374 std::cout << " Trust-region radius after update: " << del << std::endl;
375 std::cout << std::endl;
376 }
377 /***************************************************************************************************/
378 // FINISH STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
379 /***************************************************************************************************/
380 }
381
382 virtual void run( Vector<Real> &s, // Step (to be computed)
383 Real &snorm, // Step norm (to be computed)
384 int &iflag, // Exit flag (to be computed)
385 int &iter, // Iteration count (to be computed)
386 const Real del, // Trust-region radius
387 TrustRegionModel<Real> &model ) = 0; // Trust-region model
388
389 void setPredictedReduction(const Real pRed) {
390 pRed_ = pRed;
391 }
392
393 Real getPredictedReduction(void) const {
394 return pRed_;
395 }
396};
397
398}
399
401
402#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of enums for trust region algorithms.
Contains definitions of custom data types in ROL.
Provides the interface to apply upper and lower bound constraints.
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
bool isActivated(void) const
Check if bounds are on.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Provides the interface to evaluate trust-region model functions.
virtual void updatePredictedReduction(Real &pred, const Vector< Real > &s)
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
virtual const Ptr< const Vector< Real > > getGradient(void) const
virtual void updateActualReduction(Real &ared, const Vector< Real > &s)
virtual Real value(const Vector< Real > &s, Real &tol)
Compute value.
Provides interface for and implements trust-region subproblem solvers.
ETrustRegionModel TRmodel_
Real mu_
Post-Smoothing tolerance for projected methods.
Real getPredictedReduction(void) const
ROL::Ptr< Vector< Real > > xtmp_
TrustRegion(ROL::ParameterList &parlist)
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
Real alpha_init_
Initial line-search parameter for projected methods.
ROL::Ptr< Vector< Real > > dual_
ROL::Ptr< Vector< Real > > prim_
virtual void run(Vector< Real > &s, Real &snorm, int &iflag, int &iter, const Real del, TrustRegionModel< Real > &model)=0
std::vector< bool > useInexact_
void setPredictedReduction(const Real pRed)
Real beta_
Post-Smoothing rate for projected methods.
virtual void update(Vector< Real > &x, Real &fnew, Real &del, int &nfval, int &ngrad, ETrustRegionFlag &flagTR, const Vector< Real > &s, const Real snorm, const Real fold, const Vector< Real > &g, int iter, Objective< Real > &obj, BoundConstraint< Real > &bnd, TrustRegionModel< Real > &model)
int max_fval_
Maximum function evaluations in line-search for projected methods.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
Definition: ROL_Vector.hpp:226
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual Real dot(const Vector &x) const =0
Compute where .
ETrustRegionModel StringToETrustRegionModel(std::string s)
Real ROL_OVERFLOW(void)
Platform-dependent maximum double.
Definition: ROL_Types.hpp:102
@ TRUSTREGION_FLAG_POSPREDNEG
@ TRUSTREGION_FLAG_NPOSPREDNEG
@ TRUSTREGION_FLAG_QMINSUFDEC
@ TRUSTREGION_FLAG_NPOSPREDPOS
@ TRUSTREGION_FLAG_SUCCESS
@ TRUSTREGION_MODEL_KELLEYSACHS