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 
113  BundleStep(Teuchos::ParameterList &parlist)
114  : bundle_(Teuchos::null), lineSearch_(Teuchos::null),
115  QPiter_(0), QPmaxit_(0), QPtol_(0.), step_flag_(0),
116  y_(Teuchos::null), linErrNew_(0.), valueNew_(0.),
117  aggSubGradNew_(Teuchos::null), aggSubGradOldNorm_(0.),
118  aggLinErrNew_(0.), aggLinErrOld_(0.), aggDistMeasNew_(0.),
119  T_(0.), tol_(0.), m1_(0.), m2_(0.), m3_(0.), nu_(0.),
120  ls_maxit_(0), first_print_(true), isConvex_(false),
121  ftol_(ROL_EPSILON) {
122  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
123  state->searchSize = parlist.sublist("Step").sublist("Bundle").get("Initial Trust-Region Parameter", 1.e3);
124  T_ = parlist.sublist("Step").sublist("Bundle").get("Maximum Trust-Region Parameter", 1.e8);
125  tol_ = parlist.sublist("Step").sublist("Bundle").get("Epsilon Solution Tolerance", 1.e-6);
126  m1_ = parlist.sublist("Step").sublist("Bundle").get("Upper Threshold for Serious Step", 0.1);
127  m2_ = parlist.sublist("Step").sublist("Bundle").get("Lower Threshold for Serious Step", 0.2);
128  m3_ = parlist.sublist("Step").sublist("Bundle").get("Upper Threshold for Null Step", 0.9);
129  nu_ = parlist.sublist("Step").sublist("Bundle").get("Tolerance for Trust-Region Parameter", 1.e-3);
130 
131  // Initialize bundle
132  Real coeff = parlist.sublist("Step").sublist("Bundle").get("Distance Measure Coefficient", 0.0);
133  unsigned maxSize = parlist.sublist("Step").sublist("Bundle").get("Maximum Bundle Size", 200);
134  unsigned remSize = parlist.sublist("Step").sublist("Bundle").get("Removal Size for Bundle Update", 2);
135  if ( parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Solver",0) == 1 ) {
136  //bundle_ = Teuchos::rcp(new Bundle_TT<Real>(maxSize,coeff,remSize));
137  bundle_ = Teuchos::rcp(new Bundle<Real>(maxSize,coeff,remSize));
138  }
139  else {
140  bundle_ = Teuchos::rcp(new Bundle<Real>(maxSize,coeff,remSize));
141  }
142  isConvex_ = ((coeff == 0.0) ? true : false);
143 
144  // Initialize QP solver
145  QPtol_ = parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Tolerance", 1.e-8);
146  QPmaxit_ = parlist.sublist("Step").sublist("Bundle").get("Cutting Plane Iteration Limit", 1000);
147 
148  // Initialize Line Search
149  ls_maxit_
150  = parlist.sublist("Step").sublist("Line Search").get("Maximum Number of Function Evaluations",20);
151  if ( !isConvex_ ) {
152  lineSearch_ = LineSearchFactory<Real>(parlist);
153  }
154  }
155 
156  void initialize( Vector<Real> &x, const Vector<Real> &g,
158  AlgorithmState<Real> &algo_state ) {
159  // Call default initializer, but maintain current searchSize
160  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
161  Real searchSize = state->searchSize;
162  Step<Real>::initialize(x,x,g,obj,con,algo_state);
163  state->searchSize = searchSize;
164  // Initialize bundle
165  bundle_->initialize(*(state->gradientVec));
166  // Initialize storage for updated iterate
167  y_ = x.clone();
168  // Initialize storage for aggregate subgradients
169  aggSubGradNew_ = g.clone();
170  aggSubGradOldNorm_ = algo_state.gnorm;
171  // Initialize line search
172  if ( !isConvex_ ) {
173  lineSearch_->initialize(x,x,g,obj,con);
174  }
175  }
176 
178  BoundConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
179  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
180  first_print_ = false; // Print header only on first serious step
181  QPiter_ = (step_flag_ ? 0 : QPiter_); // Reset QPiter only on serious steps
182  Real v = 0.0, l = 0.0, u = T_, gd = 0.0; // Scalar storage
183  bool flag = true;
184  while (flag) {
185  /*************************************************************/
186  /******** Solve Dual Cutting Plane QP Problem ****************/
187  /*************************************************************/
188  QPiter_ += bundle_->solveDual(state->searchSize,QPmaxit_,QPtol_); // Solve QP subproblem
189  bundle_->aggregate(*aggSubGradNew_,aggLinErrNew_,aggDistMeasNew_); // Compute aggregate info
190  algo_state.aggregateGradientNorm = aggSubGradNew_->norm(); // Aggregate subgradient norm
191  /*************************************************************/
192  /******** Construct Cutting Plane Solution *******************/
193  /*************************************************************/
194  v = -state->searchSize*std::pow(algo_state.aggregateGradientNorm,2.0)-aggLinErrNew_; // CP objective value
195  s.set(aggSubGradNew_->dual()); s.scale(-state->searchSize); // CP solution
196  algo_state.snorm = state->searchSize*algo_state.aggregateGradientNorm; // Step norm
197  /*************************************************************/
198  /******** Decide Whether Step is Serious or Null *************/
199  /*************************************************************/
200  if (std::max(algo_state.aggregateGradientNorm,aggLinErrNew_) <= tol_) {
201  // Current iterate is already epsilon optimal!
202  s.zero(); algo_state.snorm = 0.0;
203  flag = false;
204  step_flag_ = 1;
205  algo_state.flag = true;
206  break;
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  bool SS1 = (valueNew_-algo_state.value < m1_*v);
221  //bool NS1 = (valueNew_-algo_state.value >= m1_*v);
222  bool NS2a = (bundle_->computeAlpha(algo_state.snorm,linErrNew_) <= m3_*aggLinErrOld_);
223  bool NS2b = (std::abs(algo_state.value-valueNew_) <= aggSubGradOldNorm_ + aggLinErrOld_);
224  if ( isConvex_ ) {
225  /************* Convex objective ****************/
226  if ( SS1 ) {
227  bool SS2 = (gd >= m2_*v || state->searchSize >= T_-nu_);
228  if ( SS2 ) { // Serious Step
229  step_flag_ = 1;
230  flag = false;
231  break;
232  }
233  else { // Increase trust-region radius
234  l = state->searchSize;
235  state->searchSize = 0.5*(u+l);
236  }
237  }
238  else {
239  if ( NS2a || NS2b ) { // Null step
240  s.zero(); algo_state.snorm = 0.0;
241  step_flag_ = 0;
242  flag = false;
243  break;
244  }
245  else { // Decrease trust-region radius
246  u = state->searchSize;
247  state->searchSize = 0.5*(u+l);
248  }
249  }
250  }
251  else {
252  /************* Nonconvex objective *************/
253  bool NS3 = (gd - bundle_->computeAlpha(algo_state.snorm,linErrNew_) >= m2_*v);
254  if ( SS1 ) { // Serious step
255  step_flag_ = 1;
256  flag = false;
257  break;
258  }
259  else {
260  if ( NS2a || NS2b ) {
261  if ( NS3 ) { // Null step
262  s.zero();
263  step_flag_ = 0;
264  flag = false;
265  break;
266  }
267  else {
268  if ( NS2b ) { // Line search
269  Real alpha = 0.0;
270  int ls_nfval = 0, ls_ngrad = 0;
271  lineSearch_->run(alpha,valueNew_,ls_nfval,ls_ngrad,gd,s,x,obj,con);
272  if ( ls_nfval == ls_maxit_ ) { // Null step
273  s.zero();
274  step_flag_ = 0;
275  flag = false;
276  break;
277  }
278  else { // Serious step
279  s.scale(alpha);
280  step_flag_ = 1;
281  flag = false;
282  break;
283  }
284  }
285  else { // Decrease trust-region radius
286  u = state->searchSize;
287  state->searchSize = 0.5*(u+l);
288  }
289  }
290  }
291  else { // Decrease trust-region radius
292  u = state->searchSize;
293  state->searchSize = 0.5*(u+l);
294  }
295  }
296  }
297  }
298  } // End While
299  /*************************************************************/
300  /******** Update Algorithm State *****************************/
301  /*************************************************************/
302  algo_state.aggregateModelError = aggLinErrNew_;
303  aggSubGradOldNorm_ = algo_state.aggregateGradientNorm;
304  aggLinErrOld_ = aggLinErrNew_;
305  } // End Compute
306 
307  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj,
308  BoundConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
309  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
310  if ( !algo_state.flag ) {
311  /*************************************************************/
312  /******** Reset Bundle If Maximum Size Reached ***************/
313  /*************************************************************/
314  bundle_->reset(*aggSubGradNew_,aggLinErrNew_,algo_state.snorm);
315  /*************************************************************/
316  /******** Update Bundle with Step Information ****************/
317  /*************************************************************/
318  if ( step_flag_ ) {
319  // Serious step was taken
320  x.plus(s); // Update iterate
321  Real valueOld = algo_state.value; // Store previous objective value
322  algo_state.value = valueNew_; // Add new objective value to state
323  bundle_->update(step_flag_,valueNew_-valueOld,algo_state.snorm,*(state->gradientVec),s);
324  }
325  else {
326  // Null step was taken
327  bundle_->update(step_flag_,linErrNew_,algo_state.snorm,*(state->gradientVec),s);
328  }
329  }
330  /*************************************************************/
331  /******** Update Algorithm State *****************************/
332  /*************************************************************/
333  algo_state.iterateVec->set(x);
334  algo_state.gnorm = (state->gradientVec)->norm();
335  if ( step_flag_ ) {
336  algo_state.iter++;
337  }
338  } // End Update
339 
340  std::string printHeader( void ) const {
341  std::stringstream hist;
342  hist << " ";
343  hist << std::setw(6) << std::left << "iter";
344  hist << std::setw(15) << std::left << "value";
345  hist << std::setw(15) << std::left << "gnorm";
346  hist << std::setw(15) << std::left << "snorm";
347  hist << std::setw(10) << std::left << "#fval";
348  hist << std::setw(10) << std::left << "#grad";
349  hist << std::setw(15) << std::left << "znorm";
350  hist << std::setw(15) << std::left << "alpha";
351  hist << std::setw(15) << std::left << "TRparam";
352  hist << std::setw(10) << std::left << "QPiter";
353  hist << "\n";
354  return hist.str();
355  }
356 
357  std::string printName( void ) const {
358  std::stringstream hist;
359  hist << "\n" << "Bundle Trust-Region Algorithm \n";
360  return hist.str();
361  }
362 
363  std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
364  Teuchos::RCP<const StepState<Real> > state = Step<Real>::getStepState();
365  std::stringstream hist;
366  hist << std::scientific << std::setprecision(6);
367  if ( algo_state.iter == 0 && first_print_ ) {
368  hist << printName();
369  if ( print_header ) {
370  hist << printHeader();
371  }
372  hist << " ";
373  hist << std::setw(6) << std::left << algo_state.iter;
374  hist << std::setw(15) << std::left << algo_state.value;
375  hist << std::setw(15) << std::left << algo_state.gnorm;
376  hist << "\n";
377  }
378  if ( step_flag_ && algo_state.iter > 0 ) {
379  if ( print_header ) {
380  hist << printHeader();
381  }
382  else {
383  hist << " ";
384  hist << std::setw(6) << std::left << algo_state.iter;
385  hist << std::setw(15) << std::left << algo_state.value;
386  hist << std::setw(15) << std::left << algo_state.gnorm;
387  hist << std::setw(15) << std::left << algo_state.snorm;
388  hist << std::setw(10) << std::left << algo_state.nfval;
389  hist << std::setw(10) << std::left << algo_state.ngrad;
390  hist << std::setw(15) << std::left << algo_state.aggregateGradientNorm;
391  hist << std::setw(15) << std::left << algo_state.aggregateModelError;
392  hist << std::setw(15) << std::left << state->searchSize;
393  hist << std::setw(10) << std::left << QPiter_;
394  hist << "\n";
395  }
396  }
397  return hist.str();
398  };
399 
400 }; // class BundleStep
401 
402 } // namespace ROL
403 
404 #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 .
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:67
Teuchos::RCP< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:72
Contains definitions of custom data types in ROL.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::string printName(void) const
Print step name.
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:77
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual Teuchos::RCP< const StepState< Real > > getStepState(void) const
Get state for step object.
Definition: ROL_Step.hpp:192
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 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.
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:87
Teuchos::RCP< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:91
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.
Teuchos::RCP< Vector< Real > > aggSubGradNew_
Teuchos::RCP< LineSearch< Real > > lineSearch_
BundleStep(Teuchos::ParameterList &parlist)
Provides the interface for and implments a bundle.
Definition: ROL_Bundle.hpp:62
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118