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"
52 #include "ROL_TrustRegionTypes.hpp"
53 #include "ROL_TrustRegionModel.hpp"
54 #include "ROL_ColemanLiModel.hpp"
55 #include "ROL_KelleySachsModel.hpp"
56 
57 namespace ROL {
58 
59 template<class Real>
60 class TrustRegion {
61 private:
62 
63  Teuchos::RCP<Vector<Real> > prim_, dual_;
64 
66 
67  Real delmax_;
68  Real eta0_, eta1_, eta2_;
70  Real pRed_;
71  Real TRsafe_, eps_;
72  Real mu0_;
73 
74  std::vector<bool> useInexact_;
75 
76  Real ftol_old_;
77 
80 
81  unsigned verbosity_;
82 
83 public:
84 
85  virtual ~TrustRegion() {}
86 
87  // Constructor
88  TrustRegion( Teuchos::ParameterList &parlist )
89  : ftol_old_(ROL_OVERFLOW<Real>()), cnt_(0), verbosity_(0) {
90  // Trust-Region Parameters
91  Teuchos::ParameterList list = parlist.sublist("Step").sublist("Trust Region");
92  TRmodel_ = StringToETrustRegionModel(list.get("Subproblem Model", "Kelley-Sachs"));
93  delmax_ = list.get("Maximum Radius", static_cast<Real>(5000.0));
94  eta0_ = list.get("Step Acceptance Threshold", static_cast<Real>(0.05));
95  eta1_ = list.get("Radius Shrinking Threshold", static_cast<Real>(0.05));
96  eta2_ = list.get("Radius Growing Threshold", static_cast<Real>(0.9));
97  gamma0_ = list.get("Radius Shrinking Rate (Negative rho)", static_cast<Real>(0.0625));
98  gamma1_ = list.get("Radius Shrinking Rate (Positive rho)", static_cast<Real>(0.25));
99  gamma2_ = list.get("Radius Growing Rate", static_cast<Real>(2.5));
100  mu0_ = list.get("Sufficient Decrease Parameter", static_cast<Real>(1.e-4));
101  TRsafe_ = list.get("Safeguard Size", static_cast<Real>(100.0));
102  eps_ = TRsafe_*ROL_EPSILON<Real>();
103  // General Inexactness Information
104  Teuchos::ParameterList &glist = parlist.sublist("General");
105  useInexact_.clear();
106  useInexact_.push_back(glist.get("Inexact Objective Function", false));
107  useInexact_.push_back(glist.get("Inexact Gradient", false));
108  useInexact_.push_back(glist.get("Inexact Hessian-Times-A-Vector", false));
109  // Inexact Function Evaluation Information
110  Teuchos::ParameterList &ilist = list.sublist("Inexact").sublist("Value");
111  scale_ = ilist.get("Tolerance Scaling", static_cast<Real>(1.e-1));
112  omega_ = ilist.get("Exponent", static_cast<Real>(0.9));
113  force_ = ilist.get("Forcing Sequence Initial Value", static_cast<Real>(1.0));
114  updateIter_ = ilist.get("Forcing Sequence Update Frequency", static_cast<int>(10));
115  forceFactor_ = ilist.get("Forcing Sequence Reduction Factor", static_cast<Real>(0.1));
116  // Get verbosity level
117  verbosity_ = glist.get("Print Verbosity", 0);
118  }
119 
120  virtual void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
121  prim_ = x.clone();
122  dual_ = g.clone();
123  }
124 
125  virtual void update( Vector<Real> &x,
126  Real &fnew,
127  Real &del,
128  int &nfval,
129  int &ngrad,
130  ETrustRegionFlag &flagTR,
131  const Vector<Real> &s,
132  const Real snorm,
133  const Real fold,
134  const Vector<Real> &g,
135  int iter,
136  Objective<Real> &obj,
138  TrustRegionModel<Real> &model ) {
139  Real tol = std::sqrt(ROL_EPSILON<Real>());
140  const Real one(1), oe4(1.e4), zero(0);
141 
142  /***************************************************************************************************/
143  // BEGIN INEXACT OBJECTIVE FUNCTION COMPUTATION
144  /***************************************************************************************************/
145  // Update inexact objective function
146  Real fold1 = fold, ftol = tol, TOL(1.e-2);
147  if ( useInexact_[0] ) {
148  if ( !(cnt_%updateIter_) && (cnt_ != 0) ) {
149  force_ *= forceFactor_;
150  }
151  Real c = scale_*std::max(TOL,std::min(one,oe4*std::max(pRed_,std::sqrt(ROL_EPSILON<Real>()))));
152  ftol = c*std::pow(std::min(eta1_,one-eta2_)
153  *std::min(std::max(pRed_,std::sqrt(ROL_EPSILON<Real>())),force_),one/omega_);
154  if ( ftol_old_ > ftol || cnt_ == 0 ) {
155  ftol_old_ = ftol;
156  fold1 = obj.value(x,ftol_old_);
157  }
158  cnt_++;
159  }
160  // Evaluate objective function at new iterate
161  prim_->set(x); prim_->plus(s);
162  obj.update(*prim_);
163  fnew = obj.value(*prim_,ftol);
164 
165  nfval = 1;
166  Real aRed = fold1 - fnew;
167  /***************************************************************************************************/
168  // FINISH INEXACT OBJECTIVE FUNCTION COMPUTATION
169  /***************************************************************************************************/
170 
171  /***************************************************************************************************/
172  // BEGIN COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
173  /***************************************************************************************************/
174  // Modify Actual and Predicted Reduction According to Model
175  model.updateActualReduction(aRed,s);
177 
178  if ( verbosity_ > 0 ) {
179  std::cout << std::endl;
180  std::cout << " Computation of actual and predicted reduction" << std::endl;
181  std::cout << " Current objective function value: " << fold1 << std::endl;
182  std::cout << " New objective function value: " << fnew << std::endl;
183  std::cout << " Actual reduction: " << aRed << std::endl;
184  std::cout << " Predicted reduction: " << pRed_ << std::endl;
185  }
186 
187  // Compute Ratio of Actual and Predicted Reduction
188  Real EPS = eps_*((one > std::abs(fold1)) ? one : std::abs(fold1));
189  Real aRed_safe = aRed + EPS, pRed_safe = pRed_ + EPS;
190  Real rho(0);
191  if (((std::abs(aRed_safe) < eps_) && (std::abs(pRed_safe) < eps_)) || aRed == pRed_) {
192  rho = one;
193  flagTR = TRUSTREGION_FLAG_SUCCESS;
194  }
195  else if ( std::isnan(aRed_safe) || std::isnan(pRed_safe) ) {
196  rho = -one;
197  flagTR = TRUSTREGION_FLAG_NAN;
198  }
199  else {
200  rho = aRed_safe/pRed_safe;
201  if (pRed_safe < zero && aRed_safe > zero) {
203  }
204  else if (aRed_safe <= zero && pRed_safe > zero) {
206  }
207  else if (aRed_safe <= zero && pRed_safe < zero) {
209  }
210  else {
211  flagTR = TRUSTREGION_FLAG_SUCCESS;
212  }
213  }
214 
215  if ( verbosity_ > 0 ) {
216  std::cout << " Safeguard: " << eps_ << std::endl;
217  std::cout << " Actual reduction with safeguard: " << aRed_safe << std::endl;
218  std::cout << " Predicted reduction with safeguard: " << pRed_safe << std::endl;
219  std::cout << " Ratio of actual and predicted reduction: " << rho << std::endl;
220  std::cout << " Trust-region flag: " << flagTR << std::endl;
221  }
222  /***************************************************************************************************/
223  // FINISH COMPUTE RATIO OF ACTUAL AND PREDICTED REDUCTION
224  /***************************************************************************************************/
225 
226 
227  /***************************************************************************************************/
228  // BEGIN CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
229  /***************************************************************************************************/
230  bool decr = true;
232  if ( rho >= eta0_ && (std::abs(aRed_safe) > eps_) ) {
233  // Compute Criticality Measure || x - P( x - g ) ||
234  prim_->set(x);
235  prim_->axpy(-one,g.dual());
236  bnd.project(*prim_);
237  prim_->scale(-one);
238  prim_->plus(x);
239  Real pgnorm = prim_->norm();
240  // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
241  prim_->set(g.dual());
242  bnd.pruneActive(*prim_,g,x);
243  Real lam = std::min(one, del/prim_->norm());
244  prim_->scale(-lam);
245  prim_->plus(x);
246  bnd.project(*prim_);
247  prim_->scale(-one);
248  prim_->plus(x);
249  pgnorm *= prim_->norm();
250  // Sufficient decrease?
251  decr = ( aRed_safe >= mu0_*eta0_*pgnorm );
252  flagTR = (!decr ? TRUSTREGION_FLAG_QMINSUFDEC : flagTR);
253 
254  if ( verbosity_ > 0 ) {
255  std::cout << " Decrease lower bound (constraints): " << 0.1*eta0_*pgnorm << std::endl;
256  std::cout << " Trust-region flag (constraints): " << flagTR << std::endl;
257  std::cout << " Is step feasible: " << bnd.isFeasible(x) << std::endl;
258  }
259  }
260  }
261  /***************************************************************************************************/
262  // FINISH CHECK SUFFICIENT DECREASE FOR BOUND CONSTRAINED PROBLEMS
263  /***************************************************************************************************/
264 
265  /***************************************************************************************************/
266  // BEGIN STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
267  /***************************************************************************************************/
268  if ( verbosity_ > 0 ) {
269  std::cout << " Norm of step: " << snorm << std::endl;
270  std::cout << " Trust-region radius before update: " << del << std::endl;
271  }
272  if ((rho < eta0_ && flagTR == TRUSTREGION_FLAG_SUCCESS) || flagTR >= 2 || !decr ) { // Step Rejected
273  fnew = fold1;
274  if (rho < zero) { // Negative reduction, interpolate to find new trust-region radius
275  Real gs(0);
276  if ( bnd.isActivated() ) {
277  model.dualTransform(*dual_, *model.getGradient());
278  gs = dual_->dot(s.dual());
279  }
280  else {
281  gs = g.dot(s.dual());
282  }
283  Real modelVal = model.value(s,tol);
284  modelVal += fold1;
285  Real theta = (one-eta2_)*gs/((one-eta2_)*(fold1+gs)+eta2_*modelVal-fnew);
286  del = std::min(gamma1_*std::min(snorm,del),std::max(gamma0_,theta)*del);
287  }
288  else { // Shrink trust-region radius
289  del = gamma1_*std::min(snorm,del);
290  }
291  obj.update(x,true,iter);
292  }
293  else if ((rho >= eta0_ && flagTR != TRUSTREGION_FLAG_NPOSPREDNEG) ||
294  (flagTR == TRUSTREGION_FLAG_POSPREDNEG)) { // Step Accepted
295  x.plus(s);
296  obj.update(x,true,iter);
297  if (rho >= eta2_) { // Increase trust-region radius
298  del = std::min(gamma2_*del,delmax_);
299  }
300  }
301  if ( verbosity_ > 0 ) {
302  std::cout << " Trust-region radius after update: " << del << std::endl;
303  std::cout << std::endl;
304  }
305  /***************************************************************************************************/
306  // FINISH STEP ACCEPTANCE AND TRUST REGION RADIUS UPDATE
307  /***************************************************************************************************/
308  }
309 
310  virtual void run( Vector<Real> &s, // Step (to be computed)
311  Real &snorm, // Step norm (to be computed)
312  int &iflag, // Exit flag (to be computed)
313  int &iter, // Iteration count (to be computed)
314  const Real del, // Trust-region radius
315  TrustRegionModel<Real> &model ) = 0; // Trust-region model
316 
317  void setPredictedReduction(const Real pRed) {
318  pRed_ = pRed;
319  }
320 
321  Real getPredictedReduction(void) const {
322  return pRed_;
323  }
324 };
325 
326 }
327 
329 
330 #endif
Provides the interface to evaluate objective functions.
TrustRegion(Teuchos::ParameterList &parlist)
virtual void plus(const Vector &x)=0
Compute , where .
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.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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:74
virtual void updatePredictedReduction(Real &pred, const Vector< Real > &s)
virtual Real dot(const Vector &x) const =0
Compute where .
virtual void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=0)
Set variables to zero if they correspond to the -active set.
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
Teuchos::RCP< Vector< Real > > prim_
bool isActivated(void)
Check if bounds are on.
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:213
ETrustRegionModel StringToETrustRegionModel(std::string s)
Real getPredictedReduction(void) const
virtual Real value(const Vector< Real > &s, Real &tol)
Compute value.
virtual const Teuchos::RCP< const Vector< Real > > getGradient(void) const
void setPredictedReduction(const Real pRed)
Real ROL_OVERFLOW(void)
Platform-dependent maximum double.
Definition: ROL_Types.hpp:151
ETrustRegionModel
Enumeration of trust-region model types.
Provides the interface to apply upper and lower bound constraints.
Teuchos::RCP< Vector< Real > > dual_
virtual void updateActualReduction(Real &ared, const Vector< Real > &s)
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
ETrustRegionFlag
Enumation of flags used by trust-region solvers.
ETrustRegionModel TRmodel_
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
Contains definitions of enums for trust region algorithms.
std::vector< bool > useInexact_
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
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)