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