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" 53 #include "ROL_BoundConstraint.hpp" 54 #include "ROL_LineSearch.hpp" 55 56 #include "ROL_ParameterList.hpp" 57 #include "ROL_Ptr.hpp" 58 59 /** @ingroup step_group 60 \class ROL::BundleStep 61 \brief Provides the interface to compute bundle trust-region steps. 62 */ 63 64 namespace ROL { 65 66 template <class Real> 67 class BundleStep : public Step<Real> { 68 private: 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 85 Real linErrNew_; 86 Real valueNew_; 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 104 int ls_maxit_; 105 106 bool first_print_; 107 bool isConvex_; 108 109 Real ftol_; 110 111 int verbosity_; 112 113 public: 114 115 using Step<Real>::initialize; 116 using Step<Real>::compute; 117 using Step<Real>::update; 118 BundleStep(ROL::ParameterList & parlist)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), 123 aggSubGradNew_(ROL::nullPtr), aggSubGradOldNorm_(0), 124 aggLinErrNew_(0), aggLinErrOld_(0), aggDistMeasNew_(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 157 ls_maxit_ 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 initialize(Vector<Real> & x,const Vector<Real> & g,Objective<Real> & obj,BoundConstraint<Real> & con,AlgorithmState<Real> & algo_state)167 void initialize( Vector<Real> &x, const Vector<Real> &g, 168 Objective<Real> &obj, BoundConstraint<Real> &con, 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 compute(Vector<Real> & s,const Vector<Real> & x,Objective<Real> & obj,BoundConstraint<Real> & con,AlgorithmState<Real> & algo_state)188 void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj, 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 /*************************************************************/ 385 algo_state.aggregateModelError = aggLinErrNew_; 386 aggSubGradOldNorm_ = algo_state.aggregateGradientNorm; 387 aggLinErrOld_ = aggLinErrNew_; 388 } // End Compute 389 update(Vector<Real> & x,const Vector<Real> & s,Objective<Real> & obj,BoundConstraint<Real> & con,AlgorithmState<Real> & algo_state)390 void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, 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 printHeader(void) const425 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 printName(void) const442 std::string printName( void ) const { 443 std::stringstream hist; 444 hist << "\n" << "Bundle Trust-Region Algorithm \n"; 445 return hist.str(); 446 } 447 print(AlgorithmState<Real> & algo_state,bool print_header=false) const448 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 490