ROL
ROL_TrustRegion.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_TRUSTREGION_H
11 #define ROL_TRUSTREGION_H
12 
17 #include "ROL_Types.hpp"
18 #include "ROL_TrustRegionTypes.hpp"
19 #include "ROL_TrustRegionModel.hpp"
20 #include "ROL_ColemanLiModel.hpp"
21 #include "ROL_KelleySachsModel.hpp"
22 
23 namespace ROL {
24 
25 template<class Real>
26 class TrustRegion {
27 private:
28 
29  ROL::Ptr<Vector<Real> > prim_, dual_, xtmp_;
30 
32 
33  Real eta0_, eta1_, eta2_;
35  Real pRed_;
36  Real TRsafe_, eps_;
37  Real mu0_;
38 
39  std::vector<bool> useInexact_;
40 
41  Real ftol_old_;
42 
45 
46  unsigned verbosity_;
47 
48  // POST SMOOTHING PARAMETERS
49  Real alpha_init_;
50  int max_fval_;
51  Real mu_;
52  Real beta_;
53 
54 public:
55 
56  virtual ~TrustRegion() {}
57 
58  // Constructor
59  TrustRegion( ROL::ParameterList &parlist )
60  : pRed_(0), ftol_old_(ROL_OVERFLOW<Real>()), cnt_(0), verbosity_(0) {
61  // Trust-Region Parameters
62  ROL::ParameterList list = parlist.sublist("Step").sublist("Trust Region");
63  TRmodel_ = StringToETrustRegionModel(list.get("Subproblem Model", "Kelley-Sachs"));
64  eta0_ = list.get("Step Acceptance Threshold", static_cast<Real>(0.05));
65  eta1_ = list.get("Radius Shrinking Threshold", static_cast<Real>(0.05));
66  eta2_ = list.get("Radius Growing Threshold", static_cast<Real>(0.9));
67  gamma0_ = list.get("Radius Shrinking Rate (Negative rho)", static_cast<Real>(0.0625));
68  gamma1_ = list.get("Radius Shrinking Rate (Positive rho)", static_cast<Real>(0.25));
69  gamma2_ = list.get("Radius Growing Rate", static_cast<Real>(2.5));
70  mu0_ = list.get("Sufficient Decrease Parameter", static_cast<Real>(1.e-4));
71  TRsafe_ = list.get("Safeguard Size", static_cast<Real>(100.0));
72  eps_ = TRsafe_*ROL_EPSILON<Real>();
73  // General Inexactness Information
74  ROL::ParameterList &glist = parlist.sublist("General");
75  useInexact_.clear();
76  useInexact_.push_back(glist.get("Inexact Objective Function", false));
77  useInexact_.push_back(glist.get("Inexact Gradient", false));
78  useInexact_.push_back(glist.get("Inexact Hessian-Times-A-Vector", false));
79  // Inexact Function Evaluation Information
80  ROL::ParameterList &ilist = list.sublist("Inexact").sublist("Value");
81  scale_ = ilist.get("Tolerance Scaling", static_cast<Real>(1.e-1));
82  omega_ = ilist.get("Exponent", static_cast<Real>(0.9));
83  force_ = ilist.get("Forcing Sequence Initial Value", static_cast<Real>(1.0));
84  updateIter_ = ilist.get("Forcing Sequence Update Frequency", static_cast<int>(10));
85  forceFactor_ = ilist.get("Forcing Sequence Reduction Factor", static_cast<Real>(0.1));
86  // Get verbosity level
87  verbosity_ = glist.get("Print Verbosity", 0);
88  // Post-smoothing parameters
89  max_fval_ = list.sublist("Post-Smoothing").get("Function Evaluation Limit", 20);
90  alpha_init_ = list.sublist("Post-Smoothing").get("Initial Step Size", static_cast<Real>(1));
91  mu_ = list.sublist("Post-Smoothing").get("Tolerance", static_cast<Real>(0.9999));
92  beta_ = list.sublist("Post-Smoothing").get("Rate", static_cast<Real>(0.01));
93  }
94 
95  virtual void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
96  prim_ = x.clone();
97  dual_ = g.clone();
98  xtmp_ = x.clone();
99  }
100 
101  virtual void update( Vector<Real> &x,
102  Real &fnew,
103  Real &del,
104  int &nfval,
105  int &ngrad,
106  ETrustRegionFlag &flagTR,
107  const Vector<Real> &s,
108  const Real snorm,
109  const Real fold,
110  const Vector<Real> &g,
111  int iter,
112  Objective<Real> &obj,
114  TrustRegionModel<Real> &model ) {
115  Real tol = std::sqrt(ROL_EPSILON<Real>());
116  const Real one(1), zero(0);
117 
118  /***************************************************************************************************/
119  // BEGIN INEXACT OBJECTIVE FUNCTION COMPUTATION
120  /***************************************************************************************************/
121  // Update inexact objective function
122  Real fold1 = fold, ftol = tol; // TOL(1.e-2);
123  if ( useInexact_[0] ) {
124  if ( !(cnt_%updateIter_) && (cnt_ != 0) ) {
125  force_ *= forceFactor_;
126  }
127  //const Real oe4(1e4);
128  //Real c = scale_*std::max(TOL,std::min(one,oe4*std::max(pRed_,std::sqrt(ROL_EPSILON<Real>()))));
129  //ftol = c*std::pow(std::min(eta1_,one-eta2_)
130  // *std::min(std::max(pRed_,std::sqrt(ROL_EPSILON<Real>())),force_),one/omega_);
131  //if ( ftol_old_ > ftol || cnt_ == 0 ) {
132  // ftol_old_ = ftol;
133  // fold1 = obj.value(x,ftol_old_);
134  //}
135  //cnt_++;
136  Real eta = static_cast<Real>(0.999)*std::min(eta1_,one-eta2_);
137  ftol = scale_*std::pow(eta*std::min(pRed_,force_),one/omega_);
138  ftol_old_ = ftol;
139  fold1 = obj.value(x,ftol_old_);
140  cnt_++;
141  }
142  // Evaluate objective function at new iterate
143  prim_->set(x); prim_->plus(s);
144  if (bnd.isActivated()) {
145  bnd.project(*prim_);
146  }
147  obj.update(*prim_);
148  fnew = obj.value(*prim_,ftol);
149 
150  nfval = 1;
151  Real aRed = fold1 - fnew;
152  /***************************************************************************************************/
153  // FINISH INEXACT OBJECTIVE FUNCTION COMPUTATION
154  /***************************************************************************************************/
155 
156  /***************************************************************************************************/
157  // BEGIN COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
158  /***************************************************************************************************/
159  // Modify Actual and Predicted Reduction According to Model
160  model.updateActualReduction(aRed,s);
162 
163  if ( verbosity_ > 0 ) {
164  std::cout << std::endl;
165  std::cout << " Computation of actual and predicted reduction" << std::endl;
166  std::cout << " Current objective function value: " << fold1 << std::endl;
167  std::cout << " New objective function value: " << fnew << std::endl;
168  std::cout << " Actual reduction: " << aRed << std::endl;
169  std::cout << " Predicted reduction: " << pRed_ << std::endl;
170  }
171 
172  // Compute Ratio of Actual and Predicted Reduction
173  Real EPS = eps_*((one > std::abs(fold1)) ? one : std::abs(fold1));
174  Real aRed_safe = aRed + EPS, pRed_safe = pRed_ + EPS;
175  Real rho(0);
176  if (((std::abs(aRed_safe) < eps_) && (std::abs(pRed_safe) < eps_)) || aRed == pRed_) {
177  rho = one;
178  flagTR = TRUSTREGION_FLAG_SUCCESS;
179  }
180  else if ( std::isnan(aRed_safe) || std::isnan(pRed_safe) ) {
181  rho = -one;
182  flagTR = TRUSTREGION_FLAG_NAN;
183  }
184  else {
185  rho = aRed_safe/pRed_safe;
186  if (pRed_safe < zero && aRed_safe > zero) {
188  }
189  else if (aRed_safe <= zero && pRed_safe > zero) {
191  }
192  else if (aRed_safe <= zero && pRed_safe < zero) {
194  }
195  else {
196  flagTR = TRUSTREGION_FLAG_SUCCESS;
197  }
198  }
199 
200  if ( verbosity_ > 0 ) {
201  std::cout << " Safeguard: " << eps_ << std::endl;
202  std::cout << " Actual reduction with safeguard: " << aRed_safe << std::endl;
203  std::cout << " Predicted reduction with safeguard: " << pRed_safe << std::endl;
204  std::cout << " Ratio of actual and predicted reduction: " << rho << std::endl;
205  std::cout << " Trust-region flag: " << flagTR << std::endl;
206  }
207  /***************************************************************************************************/
208  // FINISH COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
209  /***************************************************************************************************/
210 
211  /***************************************************************************************************/
212  // BEGIN CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
213  /***************************************************************************************************/
214  bool decr = true;
216  if ( rho >= eta0_ && (std::abs(aRed_safe) > eps_) ) {
217  // Compute Criticality Measure || x - P( x - g ) ||
218  prim_->set(x);
219  prim_->axpy(-one,g.dual());
220  bnd.project(*prim_);
221  prim_->scale(-one);
222  prim_->plus(x);
223  Real pgnorm = prim_->norm();
224  // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
225  prim_->set(g.dual());
226  bnd.pruneActive(*prim_,g,x);
227  Real lam = std::min(one, del/prim_->norm());
228  prim_->scale(-lam);
229  prim_->plus(x);
230  bnd.project(*prim_);
231  prim_->scale(-one);
232  prim_->plus(x);
233  pgnorm *= prim_->norm();
234  // Sufficient decrease?
235  decr = ( aRed_safe >= mu0_*pgnorm );
236  flagTR = (!decr ? TRUSTREGION_FLAG_QMINSUFDEC : flagTR);
237 
238  if ( verbosity_ > 0 ) {
239  std::cout << " Decrease lower bound (constraints): " << mu0_*pgnorm << std::endl;
240  std::cout << " Trust-region flag (constraints): " << flagTR << std::endl;
241  std::cout << " Is step feasible: " << bnd.isFeasible(x) << std::endl;
242  }
243  }
244  }
245  /***************************************************************************************************/
246  // FINISH CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
247  /***************************************************************************************************/
248 
249  /***************************************************************************************************/
250  // BEGIN STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
251  /***************************************************************************************************/
252  if ( verbosity_ > 0 ) {
253  std::cout << " Norm of step: " << snorm << std::endl;
254  std::cout << " Trust-region radius before update: " << del << std::endl;
255  }
256  if ((rho < eta0_ && flagTR == TRUSTREGION_FLAG_SUCCESS) || flagTR >= 2 || !decr ) { // Step Rejected
257  fnew = fold1; // This is a bug if rho < zero...
258  if (rho < zero) { // Negative reduction, interpolate to find new trust-region radius
259  Real gs(0);
260  if ( bnd.isActivated() ) {
261  model.dualTransform(*dual_, *model.getGradient());
262  gs = dual_->dot(s.dual());
263  }
264  else {
265  gs = g.dot(s.dual());
266  }
267  Real modelVal = model.value(s,tol);
268  modelVal += fold1;
269  Real theta = (one-eta2_)*gs/((one-eta2_)*(fold1+gs)+eta2_*modelVal-fnew);
270  del = std::min(gamma1_*std::min(snorm,del),std::max(gamma0_,theta)*del);
271  if ( verbosity_ > 0 ) {
272  std::cout << " Interpolation model value: " << modelVal << std::endl;
273  std::cout << " Interpolation step length: " << theta << std::endl;
274  }
275  }
276  else { // Shrink trust-region radius
277  del = gamma1_*std::min(snorm,del);
278  }
279  obj.update(x,true,iter);
280  }
281  else if ((rho >= eta0_ && flagTR != TRUSTREGION_FLAG_NPOSPREDNEG) ||
282  (flagTR == TRUSTREGION_FLAG_POSPREDNEG)) { // Step Accepted
283  // Perform line search (smoothing) to ensure decrease
285  // Compute new gradient
286  xtmp_->set(x); xtmp_->plus(s);
287  bnd.project(*xtmp_);
288  obj.update(*xtmp_);
289  obj.gradient(*dual_,*xtmp_,tol); // MUST DO SOMETHING HERE WITH TOL
290  ngrad++;
291  // Compute smoothed step
292  Real alpha(1);
293  prim_->set(*xtmp_);
294  prim_->axpy(-alpha/alpha_init_,dual_->dual());
295  bnd.project(*prim_);
296  // Compute new objective value
297  obj.update(*prim_);
298  Real ftmp = obj.value(*prim_,tol); // MUST DO SOMETHING HERE WITH TOL
299  nfval++;
300  // Perform smoothing
301  int cnt = 0;
302  alpha = alpha_init_;
303  while ( (ftmp-fnew) >= mu_*aRed ) {
304  prim_->set(*xtmp_);
305  prim_->axpy(-alpha/alpha_init_,dual_->dual());
306  bnd.project(*prim_);
307  obj.update(*prim_);
308  ftmp = obj.value(*prim_,tol); // MUST DO SOMETHING HERE WITH TOL
309  nfval++;
310  if ( cnt >= max_fval_ ) {
311  break;
312  }
313  alpha *= beta_;
314  cnt++;
315  }
316  // Store objective function and iteration information
317  if (std::isnan(ftmp)) {
318  flagTR = TRUSTREGION_FLAG_NAN;
319  del = gamma1_*std::min(snorm,del);
320  rho = static_cast<Real>(-1);
321  //x.axpy(static_cast<Real>(-1),s);
322  //obj.update(x,true,iter);
323  fnew = fold1;
324  }
325  else {
326  fnew = ftmp;
327  x.set(*prim_);
328  }
329  }
330  else {
331  x.plus(s);
332  }
333  if (rho >= eta2_) { // Increase trust-region radius
334  del = gamma2_*del;
335  }
336  obj.update(x,true,iter);
337  }
338 
339  if ( verbosity_ > 0 ) {
340  std::cout << " Trust-region radius after update: " << del << std::endl;
341  std::cout << std::endl;
342  }
343  /***************************************************************************************************/
344  // FINISH STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
345  /***************************************************************************************************/
346  }
347 
348  virtual void run( Vector<Real> &s, // Step (to be computed)
349  Real &snorm, // Step norm (to be computed)
350  int &iflag, // Exit flag (to be computed)
351  int &iter, // Iteration count (to be computed)
352  const Real del, // Trust-region radius
353  TrustRegionModel<Real> &model ) = 0; // Trust-region model
354 
355  void setPredictedReduction(const Real pRed) {
356  pRed_ = pRed;
357  }
358 
359  Real getPredictedReduction(void) const {
360  return pRed_;
361  }
362 };
363 
364 }
365 
367 
368 #endif
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
Provides the interface to evaluate objective functions.
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:192
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void plus(const Vector &x)=0
Compute , where .
bool isActivated(void) const
Check if bounds are on.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
Contains definitions of custom data types in ROL.
Provides interface for and implements trust-region subproblem solvers.
Provides the interface to evaluate trust-region model functions.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual const Ptr< const Vector< Real > > getGradient(void) const
virtual void updatePredictedReduction(Real &pred, const Vector< Real > &s)
Real alpha_init_
Initial line-search parameter for projected methods.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
ROL::Ptr< Vector< Real > > xtmp_
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Real getPredictedReduction(void) const
ETrustRegionModel StringToETrustRegionModel(std::string s)
Real mu_
Post-Smoothing tolerance for projected methods.
Real beta_
Post-Smoothing rate for projected methods.
virtual Real value(const Vector< Real > &s, Real &tol)
Compute value.
int max_fval_
Maximum function evaluations in line-search for projected methods.
void setPredictedReduction(const Real pRed)
Real ROL_OVERFLOW(void)
Platform-dependent maximum double.
Definition: ROL_Types.hpp:68
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
ETrustRegionModel
Enumeration of trust-region model types.
TrustRegion(ROL::ParameterList &parlist)
Provides the interface to apply upper and lower bound constraints.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual void updateActualReduction(Real &ared, const Vector< Real > &s)
ETrustRegionFlag
Enumation of flags used by trust-region solvers.
ETrustRegionModel TRmodel_
ROL::Ptr< Vector< Real > > dual_
Contains definitions of enums for trust region algorithms.
std::vector< bool > useInexact_
virtual void run(Vector< Real > &s, Real &snorm, int &iflag, int &iter, const Real del, TrustRegionModel< Real > &model)=0
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)
ROL::Ptr< Vector< Real > > prim_