ROL
ROL_BundleStep.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_BUNDLE_STEP_H
11 #define ROL_BUNDLE_STEP_H
12 
13 #include "ROL_Bundle_AS.hpp"
14 #include "ROL_Bundle_TT.hpp"
15 #include "ROL_Types.hpp"
16 #include "ROL_Step.hpp"
17 #include "ROL_Vector.hpp"
18 #include "ROL_Objective.hpp"
19 #include "ROL_BoundConstraint.hpp"
20 #include "ROL_LineSearch.hpp"
21 
22 #include "ROL_ParameterList.hpp"
23 #include "ROL_Ptr.hpp"
24 
30 namespace ROL {
31 
32 template <class Real>
33 class BundleStep : public Step<Real> {
34 private:
35  // Bundle
36  ROL::Ptr<Bundle<Real> > bundle_; // Bundle of subgradients and linearization errors
37  ROL::Ptr<LineSearch<Real> > lineSearch_; // Line-search object for nonconvex problems
38 
39  // Dual cutting plane solution
40  unsigned QPiter_; // Number of QP solver iterations
41  unsigned QPmaxit_; // Maximum number of QP iterations
42  Real QPtol_; // QP subproblem tolerance
43 
44  // Step flag
45  int step_flag_; // Whether serious or null step
46 
47  // Additional storage
48  ROL::Ptr<Vector<Real> > y_;
49 
50  // Updated iterate storage
51  Real linErrNew_;
52  Real valueNew_;
53 
54  // Aggregate subgradients, linearizations, and distance measures
55  ROL::Ptr<Vector<Real> > aggSubGradNew_; // New aggregate subgradient
56  Real aggSubGradOldNorm_; // Old aggregate subgradient norm
57  Real aggLinErrNew_; // New aggregate linearization error
58  Real aggLinErrOld_; // Old aggregate linearization error
59  Real aggDistMeasNew_; // New aggregate distance measure
60 
61  // Algorithmic parameters
62  Real T_;
63  Real tol_;
64  Real m1_;
65  Real m2_;
66  Real m3_;
67  Real nu_;
68 
69  // Line-search parameters
70  int ls_maxit_;
71 
73  bool isConvex_;
74 
75  Real ftol_;
76 
78 
79 public:
80 
82  using Step<Real>::compute;
83  using Step<Real>::update;
84 
85  BundleStep(ROL::ParameterList &parlist)
86  : bundle_(ROL::nullPtr), lineSearch_(ROL::nullPtr),
87  QPiter_(0), QPmaxit_(0), QPtol_(0), step_flag_(0),
88  y_(ROL::nullPtr), linErrNew_(0), valueNew_(0),
89  aggSubGradNew_(ROL::nullPtr), aggSubGradOldNorm_(0),
91  T_(0), tol_(0), m1_(0), m2_(0), m3_(0), nu_(0),
92  ls_maxit_(0), first_print_(true), isConvex_(false),
93  ftol_(ROL_EPSILON<Real>()) {
94  Real zero(0), two(2), oem3(1.e-3), oem6(1.e-6), oem8(1.e-8), p1(0.1), p2(0.2), p9(0.9), oe3(1.e3), oe8(1.e8);
95  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
96  state->searchSize = parlist.sublist("Step").sublist("Bundle").get("Initial Trust-Region Parameter", oe3);
97  T_ = parlist.sublist("Step").sublist("Bundle").get("Maximum Trust-Region Parameter", oe8);
98  tol_ = parlist.sublist("Step").sublist("Bundle").get("Epsilon Solution Tolerance", oem6);
99  m1_ = parlist.sublist("Step").sublist("Bundle").get("Upper Threshold for Serious Step", p1);
100  m2_ = parlist.sublist("Step").sublist("Bundle").get("Lower Threshold for Serious Step", p2);
101  m3_ = parlist.sublist("Step").sublist("Bundle").get("Upper Threshold for Null Step", p9);
102  nu_ = parlist.sublist("Step").sublist("Bundle").get("Tolerance for Trust-Region Parameter", oem3);
103 
104  // Initialize bundle
105  Real coeff = parlist.sublist("Step").sublist("Bundle").get("Distance Measure Coefficient", zero);
106  Real omega = parlist.sublist("Step").sublist("Bundle").get("Locality Measure Coefficient", two);
107  unsigned maxSize = parlist.sublist("Step").sublist("Bundle").get("Maximum Bundle Size", 200);
108  unsigned remSize = parlist.sublist("Step").sublist("Bundle").get("Removal Size for Bundle Update", 2);
109  if ( parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Solver",0) == 1 ) {
110  bundle_ = ROL::makePtr<Bundle_TT<Real>>(maxSize,coeff,omega,remSize);
111  //bundle_ = ROL::makePtr<Bundle_AS<Real>>(maxSize,coeff,omega,remSize);
112  }
113  else {
114  bundle_ = ROL::makePtr<Bundle_AS<Real>>(maxSize,coeff,omega,remSize);
115  }
116  isConvex_ = ((coeff == zero) ? true : false);
117 
118  // Initialize QP solver
119  QPtol_ = parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Tolerance", oem8);
120  QPmaxit_ = parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Iteration Limit", 1000);
121 
122  // Initialize Line Search
123  ls_maxit_
124  = parlist.sublist("Step").sublist("Line Search").get("Maximum Number of Function Evaluations",20);
125  if ( !isConvex_ ) {
126  lineSearch_ = LineSearchFactory<Real>(parlist);
127  }
128 
129  // Get verbosity level
130  verbosity_ = parlist.sublist("General").get("Print Verbosity", 0);
131  }
132 
133  void initialize( Vector<Real> &x, const Vector<Real> &g,
135  AlgorithmState<Real> &algo_state ) {
136  // Call default initializer, but maintain current searchSize
137  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
138  Real searchSize = state->searchSize;
139  Step<Real>::initialize(x,x,g,obj,con,algo_state);
140  state->searchSize = searchSize;
141  // Initialize bundle
142  bundle_->initialize(*(state->gradientVec));
143  // Initialize storage for updated iterate
144  y_ = x.clone();
145  // Initialize storage for aggregate subgradients
146  aggSubGradNew_ = g.clone();
147  aggSubGradOldNorm_ = algo_state.gnorm;
148  // Initialize line search
149  if ( !isConvex_ ) {
150  lineSearch_->initialize(x,x,g,obj,con);
151  }
152  }
153 
155  BoundConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
156  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
157  first_print_ = false; // Print header only on first serious step
158  QPiter_ = (step_flag_==1 ? 0 : QPiter_); // Reset QPiter only on serious steps
159  Real v(0), l(0), u = T_, gd(0); // Scalar storage
160  Real zero(0), two(2), half(0.5);
161  bool flag = true;
162  while (flag) {
163  /*************************************************************/
164  /******** Solve Dual Cutting Plane QP Problem ****************/
165  /*************************************************************/
166  QPiter_ += bundle_->solveDual(state->searchSize,QPmaxit_,QPtol_); // Solve QP subproblem
167  bundle_->aggregate(*aggSubGradNew_,aggLinErrNew_,aggDistMeasNew_); // Compute aggregate info
168  algo_state.aggregateGradientNorm = aggSubGradNew_->norm(); // Aggregate subgradient norm
169  if (verbosity_ > 0) {
170  std::cout << std::endl;
171  std::cout << " Computation of aggregrate quantities" << std::endl;
172  std::cout << " Aggregate subgradient norm: " << algo_state.aggregateGradientNorm << std::endl;
173  std::cout << " Aggregate linearization error: " << aggLinErrNew_ << std::endl;
174  std::cout << " Aggregate distance measure: " << aggDistMeasNew_ << std::endl;
175  }
176  /*************************************************************/
177  /******** Construct Cutting Plane Solution *******************/
178  /*************************************************************/
179  v = -state->searchSize*std::pow(algo_state.aggregateGradientNorm,two)-aggLinErrNew_; // CP objective value
180  s.set(aggSubGradNew_->dual()); s.scale(-state->searchSize); // CP solution
181  algo_state.snorm = state->searchSize*algo_state.aggregateGradientNorm; // Step norm
182  if (verbosity_ > 0) {
183  std::cout << std::endl;
184  std::cout << " Solve cutting plan subproblem" << std::endl;
185  std::cout << " Cutting plan objective value: " << v << std::endl;
186  std::cout << " Norm of computed step: " << algo_state.snorm << std::endl;
187  std::cout << " 'Trust-region' radius: " << state->searchSize << std::endl;
188  }
189  /*************************************************************/
190  /******** Decide Whether Step is Serious or Null *************/
191  /*************************************************************/
192  if (std::max(algo_state.aggregateGradientNorm,aggLinErrNew_) <= tol_) {
193  // Current iterate is already epsilon optimal!
194  s.zero(); algo_state.snorm = zero;
195  flag = false;
196  step_flag_ = 1;
197  algo_state.flag = true;
198  break;
199  }
200  else if (std::isnan(algo_state.aggregateGradientNorm)
201  || std::isnan(aggLinErrNew_)
202  || (std::isnan(aggDistMeasNew_) && !isConvex_)) {
203  s.zero(); algo_state.snorm = zero;
204  flag = false;
205  step_flag_ = 2;
206  algo_state.flag = true;
207  }
208  else {
209  // Current iterate is not epsilon optimal.
210  y_->set(x); y_->plus(s); // y is the candidate iterate
211  obj.update(*y_,true,algo_state.iter); // Update objective at y
212  valueNew_ = obj.value(*y_,ftol_); // Compute objective value at y
213  algo_state.nfval++;
214  obj.gradient(*(state->gradientVec),*y_,ftol_); // Compute objective (sub)gradient at y
215  algo_state.ngrad++;
216  // Compute new linearization error and distance measure
217  gd = s.dot(state->gradientVec->dual());
218  linErrNew_ = algo_state.value - (valueNew_ - gd); // Linearization error
219  // Determine whether to take a serious or null step
220  Real eps = static_cast<Real>(10)*ROL_EPSILON<Real>();
221  Real del = eps*std::max(static_cast<Real>(1),std::abs(algo_state.value));
222  Real Df = (valueNew_ - algo_state.value) - del;
223  Real Dm = v - del;
224  bool SS1 = false;
225  if (std::abs(Df) < eps && std::abs(Dm) < eps) {
226  SS1 = true;
227  }
228  else {
229  SS1 = (Df < m1_*Dm);
230  }
231  //bool SS1 = (valueNew_-algo_state.value < m1_*v);
232  //bool NS1 = (valueNew_-algo_state.value >= m1_*v);
233  bool NS2a = (bundle_->computeAlpha(algo_state.snorm,linErrNew_) <= m3_*aggLinErrOld_);
234  bool NS2b = (std::abs(algo_state.value-valueNew_) <= aggSubGradOldNorm_ + aggLinErrOld_);
235  if (verbosity_ > 0) {
236  std::cout << std::endl;
237  std::cout << " Check for serious/null step" << std::endl;
238  std::cout << " Serious step test SS(i): " << SS1 << std::endl;
239  std::cout << " -> Left hand side: " << valueNew_-algo_state.value << std::endl;
240  std::cout << " -> Right hand side: " << m1_*v << std::endl;
241  std::cout << " Null step test NS(iia): " << NS2a << std::endl;
242  std::cout << " -> Left hand side: " << bundle_->computeAlpha(algo_state.snorm,linErrNew_) << std::endl;
243  std::cout << " -> Right hand side: " << m3_*aggLinErrOld_ << std::endl;
244  std::cout << " Null step test NS(iib): " << NS2b << std::endl;
245  std::cout << " -> Left hand side: " << std::abs(algo_state.value-valueNew_) << std::endl;
246  std::cout << " -> Right hand side: " << aggSubGradOldNorm_ + aggLinErrOld_ << std::endl;
247  }
248  if ( isConvex_ ) {
249  /************* Convex objective ****************/
250  if ( SS1 ) {
251  bool SS2 = (gd >= m2_*v || state->searchSize >= T_-nu_);
252  if (verbosity_ > 0) {
253  std::cout << " Serious step test SS(iia): " << (gd >= m2_*v) << std::endl;
254  std::cout << " -> Left hand side: " << gd << std::endl;
255  std::cout << " -> Right hand side: " << m2_*v << std::endl;
256  std::cout << " Serious step test SS(iia): " << (state->searchSize >= T_-nu_) << std::endl;
257  std::cout << " -> Left hand side: " << state->searchSize << std::endl;
258  std::cout << " -> Right hand side: " << T_-nu_ << std::endl;
259  }
260  if ( SS2 ) { // Serious Step
261  step_flag_ = 1;
262  flag = false;
263  if (verbosity_ > 0) {
264  std::cout << " Serious step taken" << std::endl;
265  }
266  break;
267  }
268  else { // Increase trust-region radius
269  l = state->searchSize;
270  state->searchSize = half*(u+l);
271  if (verbosity_ > 0) {
272  std::cout << " Increase 'trust-region' radius: " << state->searchSize << std::endl;
273  }
274  }
275  }
276  else {
277  if ( NS2a || NS2b ) { // Null step
278  s.zero(); algo_state.snorm = zero;
279  step_flag_ = 0;
280  flag = false;
281  if (verbosity_ > 0) {
282  std::cout << " Null step taken" << std::endl;
283  }
284  break;
285  }
286  else { // Decrease trust-region radius
287  u = state->searchSize;
288  state->searchSize = half*(u+l);
289  if (verbosity_ > 0) {
290  std::cout << " Decrease 'trust-region' radius: " << state->searchSize << std::endl;
291  }
292  }
293  }
294  }
295  else {
296  /************* Nonconvex objective *************/
297  bool NS3 = (gd - bundle_->computeAlpha(algo_state.snorm,linErrNew_) >= m2_*v);
298  if (verbosity_ > 0) {
299  std::cout << " Null step test NS(iii): " << NS3 << std::endl;
300  std::cout << " -> Left hand side: " << gd - bundle_->computeAlpha(algo_state.snorm,linErrNew_) << std::endl;
301  std::cout << " -> Right hand side: " << m2_*v << std::endl;
302  }
303  if ( SS1 ) { // Serious step
304  step_flag_ = 1;
305  flag = false;
306  break;
307  }
308  else {
309  if ( NS2a || NS2b ) {
310  if ( NS3 ) { // Null step
311  s.zero();
312  step_flag_ = 0;
313  flag = false;
314  break;
315  }
316  else {
317  if ( NS2b ) { // Line search
318  Real alpha = zero;
319  int ls_nfval = 0, ls_ngrad = 0;
320  lineSearch_->run(alpha,valueNew_,ls_nfval,ls_ngrad,gd,s,x,obj,con);
321  if ( ls_nfval == ls_maxit_ ) { // Null step
322  s.zero();
323  step_flag_ = 0;
324  flag = false;
325  break;
326  }
327  else { // Serious step
328  s.scale(alpha);
329  step_flag_ = 1;
330  flag = false;
331  break;
332  }
333  }
334  else { // Decrease trust-region radius
335  u = state->searchSize;
336  state->searchSize = half*(u+l);
337  }
338  }
339  }
340  else { // Decrease trust-region radius
341  u = state->searchSize;
342  state->searchSize = half*(u+l);
343  }
344  }
345  }
346  }
347  } // End While
348  /*************************************************************/
349  /******** Update Algorithm State *****************************/
350  /*************************************************************/
351  algo_state.aggregateModelError = aggLinErrNew_;
354  } // End Compute
355 
356  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj,
357  BoundConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
358  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
359  state->flag = step_flag_;
360  state->SPiter = QPiter_;
361  if ( !algo_state.flag ) {
362  /*************************************************************/
363  /******** Reset Bundle If Maximum Size Reached ***************/
364  /*************************************************************/
365  bundle_->reset(*aggSubGradNew_,aggLinErrNew_,algo_state.snorm);
366  /*************************************************************/
367  /******** Update Bundle with Step Information ****************/
368  /*************************************************************/
369  if ( step_flag_==1 ) {
370  // Serious step was taken
371  x.plus(s); // Update iterate
372  Real valueOld = algo_state.value; // Store previous objective value
373  algo_state.value = valueNew_; // Add new objective value to state
374  bundle_->update(step_flag_,valueNew_-valueOld,algo_state.snorm,*(state->gradientVec),s);
375  }
376  else if ( step_flag_==0 ) {
377  // Null step was taken
378  bundle_->update(step_flag_,linErrNew_,algo_state.snorm,*(state->gradientVec),s);
379  }
380  }
381  /*************************************************************/
382  /******** Update Algorithm State *****************************/
383  /*************************************************************/
384  algo_state.iterateVec->set(x);
385  algo_state.gnorm = (state->gradientVec)->norm();
386  if ( step_flag_==1 ) {
387  algo_state.iter++;
388  }
389  } // End Update
390 
391  std::string printHeader( void ) const {
392  std::stringstream hist;
393  hist << " ";
394  hist << std::setw(6) << std::left << "iter";
395  hist << std::setw(15) << std::left << "value";
396  hist << std::setw(15) << std::left << "gnorm";
397  hist << std::setw(15) << std::left << "snorm";
398  hist << std::setw(10) << std::left << "#fval";
399  hist << std::setw(10) << std::left << "#grad";
400  hist << std::setw(15) << std::left << "znorm";
401  hist << std::setw(15) << std::left << "alpha";
402  hist << std::setw(15) << std::left << "TRparam";
403  hist << std::setw(10) << std::left << "QPiter";
404  hist << "\n";
405  return hist.str();
406  }
407 
408  std::string printName( void ) const {
409  std::stringstream hist;
410  hist << "\n" << "Bundle Trust-Region Algorithm \n";
411  return hist.str();
412  }
413 
414  std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
415  const ROL::Ptr<const StepState<Real> > state = Step<Real>::getStepState();
416  std::stringstream hist;
417  hist << std::scientific << std::setprecision(6);
418  if ( algo_state.iter == 0 && first_print_ ) {
419  hist << printName();
420  if ( print_header ) {
421  hist << printHeader();
422  }
423  hist << " ";
424  hist << std::setw(6) << std::left << algo_state.iter;
425  hist << std::setw(15) << std::left << algo_state.value;
426  hist << std::setw(15) << std::left << algo_state.gnorm;
427  hist << "\n";
428  }
429  if ( step_flag_==1 && algo_state.iter > 0 ) {
430  if ( print_header ) {
431  hist << printHeader();
432  }
433  else {
434  hist << " ";
435  hist << std::setw(6) << std::left << algo_state.iter;
436  hist << std::setw(15) << std::left << algo_state.value;
437  hist << std::setw(15) << std::left << algo_state.gnorm;
438  hist << std::setw(15) << std::left << algo_state.snorm;
439  hist << std::setw(10) << std::left << algo_state.nfval;
440  hist << std::setw(10) << std::left << algo_state.ngrad;
441  hist << std::setw(15) << std::left << algo_state.aggregateGradientNorm;
442  hist << std::setw(15) << std::left << algo_state.aggregateModelError;
443  hist << std::setw(15) << std::left << state->searchSize;
444  hist << std::setw(10) << std::left << QPiter_;
445  hist << "\n";
446  }
447  }
448  return hist.str();
449  };
450 
451 }; // class BundleStep
452 
453 } // namespace ROL
454 
455 #endif
Provides the interface to evaluate objective functions.
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Ptr< Vector< Real > > y_
virtual void plus(const Vector &x)=0
Compute , where .
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:34
Contains definitions of custom data types in ROL.
BundleStep(ROL::ParameterList &parlist)
std::string printName(void) const
Print step name.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:109
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:39
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with bound constraint.
ROL::Ptr< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:123
std::string printHeader(void) const
Print iterate header.
Provides the interface to apply upper and lower bound constraints.
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful.
std::string print(AlgorithmState< Real > &algo_state, bool print_header=false) const
Print iterate status.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step.
Provides the interface to compute bundle trust-region steps.
ROL::Ptr< Vector< Real > > aggSubGradNew_
virtual void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with bound constraint.
Definition: ROL_Step.hpp:54
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:57
ROL::Ptr< Bundle< Real > > bundle_
ROL::Ptr< LineSearch< Real > > lineSearch_
const ROL::Ptr< const StepState< Real > > getStepState(void) const
Get state for step object.
Definition: ROL_Step.hpp:177