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_HelperFunctions.hpp"
53 
54 namespace ROL {
55 
56 template<class Real>
57 class TrustRegion {
58 private:
59 
60  Teuchos::RCP<Vector<Real> > xupdate_;
61  Teuchos::RCP<Vector<Real> > Hs_;
62 
63  Real delmax_;
64  Real eta0_;
65  Real eta1_;
66  Real eta2_;
67  Real gamma0_;
68  Real gamma1_;
69  Real gamma2_;
70 
71  Real pRed_;
72 
73  Real TRsafe_;
74  Real eps_;
75 
76  std::vector<bool> useInexact_;
77 
78  Real ftol_old_;
79 
80  Real scale_;
81  Real omega_;
82  Real force_;
85  int cnt_;
86 
87  bool softUp_;
88 
89  unsigned verbosity_;
90 
91  void updateObj( Vector<Real> &x, int iter, ProjectedObjective<Real> &pObj ) {
92  if ( !softUp_ ) {
93  pObj.update(x,true,iter);
94  }
95  else {
96  pObj.update(x);
97  }
98  }
99 
100 
101 public:
102 
103  virtual ~TrustRegion() {}
104 
105  // Constructor
106  TrustRegion( Teuchos::ParameterList & parlist )
107  : ftol_old_(ROL_OVERFLOW), cnt_(0), verbosity_(0) {
108  // Unravel Parameter List
109  // Trust-Region Parameters
110  Teuchos::ParameterList list = parlist.sublist("Step").sublist("Trust Region");
111  delmax_ = list.get("Maximum Radius",5000.0);
112  eta0_ = list.get("Step Acceptance Threshold",0.05);
113  eta1_ = list.get("Radius Shrinking Threshold",0.05);
114  eta2_ = list.get("Radius Growing Threshold",0.9);
115  gamma0_ = list.get("Radius Shrinking Rate (Negative rho)",0.0625);
116  gamma1_ = list.get("Radius Shrinking Rate (Positive rho)",0.25);
117  gamma2_ = list.get("Radius Growing Rate",2.5);
118  TRsafe_ = list.get("Safeguard Size",100.0);
119  eps_ = TRsafe_*ROL_EPSILON;
120 
121  // Inexactness Information
122  useInexact_.clear();
123  useInexact_.push_back(parlist.sublist("General").get("Inexact Objective Function", false));
124  useInexact_.push_back(parlist.sublist("General").get("Inexact Gradient", false));
125  useInexact_.push_back(parlist.sublist("General").get("Inexact Hessian-Times-A-Vector", false));
126  scale_ = list.sublist("Inexact").sublist("Value").get("Tolerance Scaling",1.e-1);
127  omega_ = list.sublist("Inexact").sublist("Value").get("Exponent",0.9);
128  force_ = list.sublist("Inexact").sublist("Value").get("Forcing Sequence Initial Value",1.0);
129  updateIter_ = list.sublist("Inexact").sublist("Value").get("Forcing Sequence Update Frequency",10);
130  forceFactor_ = list.sublist("Inexact").sublist("Value").get("Forcing Sequence Reduction Factor",0.1);
131 
132  // Changing Objective Functions
133  softUp_ = parlist.sublist("General").get("Variable Objective Function",false);
134  }
135 
136  virtual void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
137  xupdate_ = x.clone();
138  Hs_ = g.clone();
139  }
140 
141  virtual void update( Vector<Real> &x,
142  Real &fnew,
143  Real &del,
144  int &nfval,
145  int &ngrad,
146  int &flagTR,
147  const Vector<Real> &s,
148  const Real snorm,
149  const Real fold,
150  const Vector<Real> &g,
151  int iter,
152  ProjectedObjective<Real> &pObj ) {
153  Real tol = std::sqrt(ROL_EPSILON);
154 
155  // Compute updated iterate vector
156  xupdate_->set(x);
157  xupdate_->axpy(1.0,s);
158  /***************************************************************************************************/
159  // BEGIN OBJECTIVE FUNCTION COMPUTATION
160  /***************************************************************************************************/
161  // Update inexact objective function
162  Real fold1 = fold, ftol = tol;
163  if ( useInexact_[0] ) {
164  if ( !(cnt_%updateIter_) && (cnt_ != 0) ) {
165  force_ *= forceFactor_;
166  }
167  Real c = scale_*std::max(1.e-2,std::min(1.0,1.e4*std::max(pRed_,std::sqrt(ROL_EPSILON))));
168  ftol = c*std::pow(std::min(eta1_,1.0-eta2_)
169  *std::min(std::max(pRed_,std::sqrt(ROL_EPSILON)),force_),1.0/omega_);
170  if ( ftol_old_ > ftol || cnt_ == 0 ) {
171  ftol_old_ = ftol;
172  fold1 = pObj.value(x,ftol_old_);
173  }
174  cnt_++;
175  }
176  // Evaluate objective function at new iterat
177  //updateObj(*xupdate_,iter,pObj);
178  //pObj.update(*xupdate_,true,iter);
179  pObj.update(*xupdate_);
180  fnew = pObj.value(*xupdate_,ftol);
181  pObj.update(x);
182 
183  nfval = 1;
184  Real aRed = fold1 - fnew;
185  /***************************************************************************************************/
186  // FINISH OBJECTIVE FUNCTION COMPUTATION
187  /***************************************************************************************************/
188  // If constraints are turned on, then compute a different predicted reduction
189  if (pObj.isConActivated()) {
190  xupdate_->set(x);
191  xupdate_->axpy(-1.0,g.dual());
192  pObj.project(*xupdate_);
193  xupdate_->axpy(-1.0,x);
194  xupdate_->scale(-1.0);
195 
196  pObj.reducedHessVec(*Hs_,s,x,xupdate_->dual(),x,tol);
197  pRed_ = -0.5*s.dot(Hs_->dual());
198 
199  Hs_->set(g);
200  pObj.pruneActive(*Hs_,xupdate_->dual(),x);
201  pRed_ -= s.dot(Hs_->dual());
202  }
203 
204  if ( verbosity_ > 0 ) {
205  std::cout << std::endl;
206  std::cout << " Computation of actual and predicted reduction" << std::endl;
207  std::cout << " Current objective function value: " << fold1 << std::endl;
208  std::cout << " New objective function value: " << fnew << std::endl;
209  std::cout << " Actual reduction: " << aRed << std::endl;
210  std::cout << " Predicted reduction: " << pRed_ << std::endl;
211  }
212 
213  // Compute Ratio of Actual and Predicted Reduction
214  aRed -= eps_*((1.0 > std::abs(fold1)) ? 1.0 : std::abs(fold1));
215  pRed_ -= eps_*((1.0 > std::abs(fold1)) ? 1.0 : std::abs(fold1));
216  Real rho = 0.0;
217  if ((std::abs(aRed) < eps_) && (std::abs(pRed_) < eps_)) {
218  rho = 1.0;
219  flagTR = 0;
220  }
221  else if ( std::isnan(aRed) || std::isnan(pRed_) ) {
222  rho = -1.0;
223  flagTR = 5;
224  }
225  else {
226  rho = aRed/pRed_;
227  if (pRed_ < 0 && aRed > 0) {
228  flagTR = 1;
229  }
230  else if (aRed <= 0 && pRed_ > 0) {
231  flagTR = 2;
232  }
233  else if (aRed <= 0 && pRed_ < 0) {
234  flagTR = 3;
235  }
236  else {
237  flagTR = 0;
238  }
239  }
240 
241  if ( verbosity_ > 0 ) {
242  std::cout << " Actual reduction with safeguard: " << aRed << std::endl;
243  std::cout << " Predicted reduction with safeguard: " << pRed_ << std::endl;
244  std::cout << " Ratio of actual and predicted reduction: " << rho << std::endl;
245  std::cout << " Trust-region flag: " << flagTR << std::endl;
246  }
247 
248  // Check Sufficient Decrease in the Reduced Quadratic Model
249  bool decr = true;
250  if ( pObj.isConActivated() && (std::abs(aRed) > eps_) ) {
251  // Compute Criticality Measure || x - P( x - g ) ||
252  xupdate_->set(x);
253  xupdate_->axpy(-1.0,g.dual());
254  pObj.project(*xupdate_);
255  xupdate_->scale(-1.0);
256  xupdate_->plus(x);
257  Real pgnorm = xupdate_->norm();
258  // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
259  xupdate_->set(g.dual());
260  pObj.pruneActive(*xupdate_,g,x);
261  Real lam = std::min(1.0, del/xupdate_->norm());
262  xupdate_->scale(-lam);
263  xupdate_->plus(x);
264  pObj.project(*xupdate_);
265  xupdate_->scale(-1.0);
266  xupdate_->plus(x);
267  pgnorm *= xupdate_->norm();
268  // Sufficient decrease?
269  decr = ( aRed >= 0.1*eta0_*pgnorm );
270  flagTR = (!decr ? 4 : flagTR);
271 
272  if ( verbosity_ > 0 ) {
273  std::cout << " Decrease lower bound (constraints): " << 0.1*eta0_*pgnorm << std::endl;
274  std::cout << " Trust-region flag (constraints): " << flagTR << std::endl;
275  std::cout << " Is step feasible: " << pObj.isFeasible(x) << std::endl;
276  }
277  }
278 
279  if ( verbosity_ > 0 ) {
280  std::cout << std::endl;
281  }
282 
283  // Accept or Reject Step and Update Trust Region
284  if ((rho < eta0_ && flagTR == 0) || flagTR >= 2 || !decr ) { // Step Rejected
285  //updateObj(x,iter,pObj);
286  //pObj.update(x,true,iter);
287  fnew = fold1;
288  if (rho < 0.0) { // Negative reduction, interpolate to find new trust-region radius
289  Real gs = s.dot(g.dual());
290  pObj.hessVec(*Hs_,s,x,tol);
291  Real modelVal = s.dot(Hs_->dual());
292  modelVal *= 0.5;
293  modelVal += gs + fold1;
294  Real theta = (1.0-eta2_)*gs/((1.0-eta2_)*(fold1+gs)+eta2_*modelVal-fnew);
295  del = std::min(gamma1_*snorm,std::max(gamma0_,theta)*del);
296  }
297  else { // Shrink trust-region radius
298  del = gamma1_*snorm;
299  }
300  }
301  else if ((rho >= eta0_ && flagTR != 3) || flagTR == 1) { // Step Accepted
302  x.axpy(1.0,s);
303  pObj.update(x,true,iter);
304  if (rho >= eta2_) { // Increase trust-region radius
305  del = std::min(gamma2_*del,delmax_);
306  }
307  }
308  else { // step rejected
309  if(softUp_) { // Variable Objective Function
310  pObj.update(x,true,iter);
311  }
312  }
313  }
314 
315  void setPredictedReduction(const Real pRed) {
316  pRed_ = pRed;
317  }
318  Real getPredictedReduction(void) const {
319  return pRed_;
320  }
321 
322  virtual void run( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
323  const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) = 0;
324 #if 0
325  // If constraints are active, then determine a feasible step
326  if ( pObj.isConActivated() && etr_ != TRUSTREGION_CAUCHYPOINT ) {
327  // Compute projected step stmp = P( x + alpha*s ) - x and xtmp = x + stmp
328  Real alpha = 1.0;
329  Teuchos::RCP<Vector<Real> > stmp = x.clone();
330  stmp->set(s);
331  stmp->scale(alpha);
332  pObj.computeProjectedStep(*stmp,x);
333  Teuchos::RCP<Vector<Real> > xtmp = x.clone();
334  xtmp->set(x);
335  xtmp->axpy(1.0,*stmp);
336  // Compute model components for alpha = 1.0
337  Real tol = std::sqrt(ROL_EPSILON);
338  Teuchos::RCP<Vector<Real> > Bs = x.clone();
339  pObj.hessVec(*Bs,*stmp,x,tol);
340  Real sBs = Bs->dot(*stmp);
341  Real gs = grad.dot(*stmp);
342  Real val = gs + 0.5*sBs;
343  Real val0 = val;
344  // Backtrack alpha until x+alpha*s is feasible
345  int cnt = 0;
346  int maxit = 10;
347  while ( val > val0 || !pObj.isFeasible(*xtmp) ) {
348  // Backtrack alpha
349  alpha *= 0.5;
350  // Computed projected step P( x + alpha*s ) - x and xtmp = x + stmp
351  stmp->set(s);
352  stmp->scale(alpha);
353  pObj.computeProjectedStep(*stmp,x);
354  xtmp->set(x);
355  xtmp->axpy(1.0,*stmp);
356  // Evaluate Model
357  val0 = val;
358  pObj.hessVec(*Bs,*stmp,x,tol);
359  sBs = Bs->dot(*stmp);
360  gs = grad.dot(*stmp);
361  val = gs + 0.5*sBs;
362  // Update iteration count
363  cnt++;
364  if ( cnt >= maxit ) { break; }
365  }
366  s.set(*stmp);
367  pRed_ = -val;
368  }
369 #endif
370 };
371 
372 }
373 
375 
376 #endif
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
Real value(const Vector< Real > &x, Real &tol)
Compute value.
TrustRegion(Teuchos::ParameterList &parlist)
Teuchos::RCP< Vector< Real > > xupdate_
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
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.
bool isFeasible(const Vector< Real > &v)
void pruneActive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
Contains definitions for helper functions in ROL.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual Real dot(const Vector &x) const =0
Compute where .
virtual void run(Vector< Real > &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector< Real > &x, const Vector< Real > &grad, const Real &gnorm, ProjectedObjective< Real > &pObj)=0
Real getPredictedReduction(void) const
Teuchos::RCP< Vector< Real > > Hs_
void updateObj(Vector< Real > &x, int iter, ProjectedObjective< Real > &pObj)
virtual void update(Vector< Real > &x, Real &fnew, Real &del, int &nfval, int &ngrad, int &flagTR, const Vector< Real > &s, const Real snorm, const Real fold, const Vector< Real > &g, int iter, ProjectedObjective< Real > &pObj)
void setPredictedReduction(const Real pRed)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void computeProjectedStep(Vector< Real > &v, const Vector< Real > &x)
void project(Vector< Real > &x)
static const double ROL_OVERFLOW
Platform-dependent maximum double.
Definition: ROL_Types.hpp:126
std::vector< bool > useInexact_
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118