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