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.
virtual void set(const Vector &x)
Set where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
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.
@ TRUSTREGION_FLAG_POSPREDNEG
@ TRUSTREGION_FLAG_NPOSPREDNEG
@ TRUSTREGION_FLAG_QMINSUFDEC
@ TRUSTREGION_FLAG_NPOSPREDPOS
@ TRUSTREGION_MODEL_KELLEYSACHS