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_COMPOSITESTEP_H 45 #define ROL_COMPOSITESTEP_H 46 47 #include "ROL_Types.hpp" 48 #include "ROL_Step.hpp" 49 #include "ROL_LAPACK.hpp" 50 #include "ROL_LinearAlgebra.hpp" 51 #include <sstream> 52 #include <iomanip> 53 54 /** \class ROL::CompositeStep 55 \brief Implements the computation of optimization steps 56 with composite-step trust-region methods. 57 */ 58 59 60 namespace ROL { 61 62 template <class Real> 63 class CompositeStep : public Step<Real> { 64 private: 65 66 // Vectors used for cloning. 67 ROL::Ptr<Vector<Real> > xvec_; 68 ROL::Ptr<Vector<Real> > gvec_; 69 ROL::Ptr<Vector<Real> > cvec_; 70 ROL::Ptr<Vector<Real> > lvec_; 71 72 // Diagnostic return flags for subalgorithms. 73 int flagCG_; 74 int flagAC_; 75 int iterCG_; 76 77 // Stopping conditions. 78 int maxiterCG_; 79 int maxiterOSS_; 80 Real tolCG_; 81 Real tolOSS_; 82 bool tolOSSfixed_; 83 84 // Tolerances and stopping conditions for subalgorithms. 85 Real lmhtol_; 86 Real qntol_; 87 Real pgtol_; 88 Real projtol_; 89 Real tangtol_; 90 Real tntmax_; 91 92 // Trust-region parameters. 93 Real zeta_; 94 Real Delta_; 95 Real penalty_; 96 Real eta_; 97 bool useConHess_; 98 99 Real ared_; 100 Real pred_; 101 Real snorm_; 102 Real nnorm_; 103 Real tnorm_; 104 105 // Output flags. 106 bool infoQN_; 107 bool infoLM_; 108 bool infoTS_; 109 bool infoAC_; 110 bool infoLS_; 111 bool infoALL_; 112 113 // Performance summary. 114 int totalIterCG_; 115 int totalProj_; 116 int totalNegCurv_; 117 int totalRef_; 118 int totalCallLS_; 119 int totalIterLS_; 120 sgn(T val)121 template <typename T> int sgn(T val) { 122 return (T(0) < val) - (val < T(0)); 123 } 124 printInfoLS(const std::vector<Real> & res) const125 void printInfoLS(const std::vector<Real> &res) const { 126 if (infoLS_) { 127 std::stringstream hist; 128 hist << std::scientific << std::setprecision(8); 129 hist << "\n Augmented System Solver:\n"; 130 hist << " True Residual\n"; 131 for (unsigned j=0; j<res.size(); j++) { 132 hist << " " << std::left << std::setw(14) << res[j] << "\n"; 133 } 134 hist << "\n"; 135 std::cout << hist.str(); 136 } 137 } 138 setTolOSS(const Real intol) const139 Real setTolOSS(const Real intol) const { 140 return tolOSSfixed_ ? tolOSS_ : intol; 141 } 142 143 public: 144 using Step<Real>::initialize; 145 using Step<Real>::compute; 146 using Step<Real>::update; 147 ~CompositeStep()148 virtual ~CompositeStep() {} 149 CompositeStep(ROL::ParameterList & parlist)150 CompositeStep( ROL::ParameterList & parlist ) : Step<Real>() { 151 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState(); 152 flagCG_ = 0; 153 flagAC_ = 0; 154 iterCG_ = 0; 155 156 ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Composite Step"); 157 158 //maxiterOSS_ = steplist.sublist("Optimality System Solver").get("Iteration Limit", 50); 159 tolOSS_ = steplist.sublist("Optimality System Solver").get("Nominal Relative Tolerance", 1e-8); 160 tolOSSfixed_ = steplist.sublist("Optimality System Solver").get("Fix Tolerance", true); 161 162 maxiterCG_ = steplist.sublist("Tangential Subproblem Solver").get("Iteration Limit", 20); 163 tolCG_ = steplist.sublist("Tangential Subproblem Solver").get("Relative Tolerance", 1e-2); 164 Delta_ = steplist.get("Initial Radius", 1e2); 165 useConHess_ = steplist.get("Use Constraint Hessian", true); 166 167 int outLvl = steplist.get("Output Level", 0); 168 169 lmhtol_ = tolOSS_; 170 qntol_ = tolOSS_; 171 pgtol_ = tolOSS_; 172 projtol_ = tolOSS_; 173 tangtol_ = tolOSS_; 174 tntmax_ = 2.0; 175 176 zeta_ = 0.8; 177 penalty_ = 1.0; 178 eta_ = 1e-8; 179 180 snorm_ = 0.0; 181 nnorm_ = 0.0; 182 tnorm_ = 0.0; 183 184 infoALL_ = false; 185 if (outLvl > 0) { 186 infoALL_ = true; 187 } 188 infoQN_ = false; 189 infoLM_ = false; 190 infoTS_ = false; 191 infoAC_ = false; 192 infoLS_ = false; 193 infoQN_ = infoQN_ || infoALL_; 194 infoLM_ = infoLM_ || infoALL_; 195 infoTS_ = infoTS_ || infoALL_; 196 infoAC_ = infoAC_ || infoALL_; 197 infoLS_ = infoLS_ || infoALL_; 198 199 totalIterCG_ = 0; 200 totalProj_ = 0; 201 totalNegCurv_ = 0; 202 totalRef_ = 0; 203 totalCallLS_ = 0; 204 totalIterLS_ = 0; 205 } 206 207 /** \brief Initialize step. 208 */ initialize(Vector<Real> & x,const Vector<Real> & g,Vector<Real> & l,const Vector<Real> & c,Objective<Real> & obj,Constraint<Real> & con,AlgorithmState<Real> & algo_state)209 void initialize( Vector<Real> &x, const Vector<Real> &g, Vector<Real> &l, const Vector<Real> &c, 210 Objective<Real> &obj, Constraint<Real> &con, 211 AlgorithmState<Real> &algo_state ) { 212 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState(); 213 ROL::Ptr<StepState<Real> > state = Step<Real>::getState(); 214 state->descentVec = x.clone(); 215 state->gradientVec = g.clone(); 216 state->constraintVec = c.clone(); 217 218 xvec_ = x.clone(); 219 gvec_ = g.clone(); 220 lvec_ = l.clone(); 221 cvec_ = c.clone(); 222 223 ROL::Ptr<Vector<Real> > ajl = gvec_->clone(); 224 ROL::Ptr<Vector<Real> > gl = gvec_->clone(); 225 226 algo_state.nfval = 0; 227 algo_state.ncval = 0; 228 algo_state.ngrad = 0; 229 230 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); 231 232 // Update objective and constraint. 233 obj.update(x,true,algo_state.iter); 234 algo_state.value = obj.value(x, zerotol); 235 algo_state.nfval++; 236 con.update(x,true,algo_state.iter); 237 con.value(*cvec_, x, zerotol); 238 algo_state.cnorm = cvec_->norm(); 239 algo_state.ncval++; 240 obj.gradient(*gvec_, x, zerotol); 241 242 // Compute gradient of Lagrangian at new multiplier guess. 243 computeLagrangeMultiplier(l, x, *gvec_, con); 244 con.applyAdjointJacobian(*ajl, l, x, zerotol); 245 gl->set(*gvec_); gl->plus(*ajl); 246 algo_state.ngrad++; 247 algo_state.gnorm = gl->norm(); 248 } 249 250 /** \brief Compute step. 251 */ compute(Vector<Real> & s,const Vector<Real> & x,const Vector<Real> & l,Objective<Real> & obj,Constraint<Real> & con,AlgorithmState<Real> & algo_state)252 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l, 253 Objective<Real> &obj, Constraint<Real> &con, 254 AlgorithmState<Real> &algo_state ) { 255 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState(); 256 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); 257 Real f = 0.0; 258 ROL::Ptr<Vector<Real> > n = xvec_->clone(); 259 ROL::Ptr<Vector<Real> > c = cvec_->clone(); 260 ROL::Ptr<Vector<Real> > t = xvec_->clone(); 261 ROL::Ptr<Vector<Real> > tCP = xvec_->clone(); 262 ROL::Ptr<Vector<Real> > g = gvec_->clone(); 263 ROL::Ptr<Vector<Real> > gf = gvec_->clone(); 264 ROL::Ptr<Vector<Real> > Wg = xvec_->clone(); 265 ROL::Ptr<Vector<Real> > ajl = gvec_->clone(); 266 267 Real f_new = 0.0; 268 ROL::Ptr<Vector<Real> > l_new = lvec_->clone(); 269 ROL::Ptr<Vector<Real> > c_new = cvec_->clone(); 270 ROL::Ptr<Vector<Real> > g_new = gvec_->clone(); 271 ROL::Ptr<Vector<Real> > gf_new = gvec_->clone(); 272 273 // Evaluate objective ... should have been stored. 274 f = obj.value(x, zerotol); 275 algo_state.nfval++; 276 // Compute gradient of objective ... should have been stored. 277 obj.gradient(*gf, x, zerotol); 278 // Evaluate constraint ... should have been stored. 279 con.value(*c, x, zerotol); 280 281 // Compute quasi-normal step. 282 computeQuasinormalStep(*n, *c, x, zeta_*Delta_, con); 283 284 // Compute gradient of Lagrangian ... should have been stored. 285 con.applyAdjointJacobian(*ajl, l, x, zerotol); 286 g->set(*gf); 287 g->plus(*ajl); 288 algo_state.ngrad++; 289 290 // Solve tangential subproblem. 291 solveTangentialSubproblem(*t, *tCP, *Wg, x, *g, *n, l, Delta_, obj, con); 292 totalIterCG_ += iterCG_; 293 294 // Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence. 295 accept(s, *n, *t, f_new, *c_new, *gf_new, *l_new, *g_new, x, l, f, *gf, *c, *g, *tCP, *Wg, obj, con, algo_state); 296 } 297 298 /** \brief Update step, if successful. 299 */ update(Vector<Real> & x,Vector<Real> & l,const Vector<Real> & s,Objective<Real> & obj,Constraint<Real> & con,AlgorithmState<Real> & algo_state)300 void update( Vector<Real> &x, Vector<Real> &l, const Vector<Real> &s, 301 Objective<Real> &obj, Constraint<Real> &con, 302 AlgorithmState<Real> &algo_state ) { 303 //ROL::Ptr<StepState<Real> > state = Step<Real>::getState(); 304 305 Real zero(0); 306 Real one(1); 307 Real two(2); 308 Real seven(7); 309 Real half(0.5); 310 Real zp9(0.9); 311 Real zp8(0.8); 312 Real em12(1e-12); 313 Real zerotol = std::sqrt(ROL_EPSILON<Real>());//zero; 314 Real ratio(zero); 315 316 ROL::Ptr<Vector<Real> > g = gvec_->clone(); 317 ROL::Ptr<Vector<Real> > ajl = gvec_->clone(); 318 ROL::Ptr<Vector<Real> > gl = gvec_->clone(); 319 ROL::Ptr<Vector<Real> > c = cvec_->clone(); 320 321 // Determine if the step gives sufficient reduction in the merit function, 322 // update the trust-region radius. 323 ratio = ared_/pred_; 324 if ((std::abs(ared_) < em12) && std::abs(pred_) < em12) { 325 ratio = one; 326 } 327 if (ratio >= eta_) { 328 x.plus(s); 329 if (ratio >= zp9) { 330 Delta_ = std::max(seven*snorm_, Delta_); 331 } 332 else if (ratio >= zp8) { 333 Delta_ = std::max(two*snorm_, Delta_); 334 } 335 obj.update(x,true,algo_state.iter); 336 con.update(x,true,algo_state.iter); 337 flagAC_ = 1; 338 } 339 else { 340 Delta_ = half*std::max(nnorm_, tnorm_); 341 obj.update(x,false,algo_state.iter); 342 con.update(x,false,algo_state.iter); 343 flagAC_ = 0; 344 } // if (ratio >= eta) 345 346 Real val = obj.value(x, zerotol); 347 algo_state.nfval++; 348 obj.gradient(*g, x, zerotol); 349 computeLagrangeMultiplier(l, x, *g, con); 350 con.applyAdjointJacobian(*ajl, l, x, zerotol); 351 gl->set(*g); gl->plus(*ajl); 352 algo_state.ngrad++; 353 con.value(*c, x, zerotol); 354 355 ROL::Ptr<StepState<Real> > state = Step<Real>::getState(); 356 state->gradientVec->set(*gl); 357 state->constraintVec->set(*c); 358 359 algo_state.value = val; 360 algo_state.gnorm = gl->norm(); 361 algo_state.cnorm = c->norm(); 362 algo_state.iter++; 363 algo_state.snorm = snorm_; 364 365 // Update algorithm state 366 //(algo_state.iterateVec)->set(x); 367 } 368 369 /** \brief Compute step for bound constraints; here only to satisfy the 370 interface requirements, does nothing, needs refactoring. 371 */ compute(Vector<Real> & s,const Vector<Real> & x,Objective<Real> & obj,BoundConstraint<Real> & con,AlgorithmState<Real> & algo_state)372 void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj, 373 BoundConstraint<Real> &con, 374 AlgorithmState<Real> &algo_state ) {} 375 376 /** \brief Update step, for bound constraints; here only to satisfy the 377 interface requirements, does nothing, needs refactoring. 378 */ update(Vector<Real> & x,const Vector<Real> & s,Objective<Real> & obj,BoundConstraint<Real> & con,AlgorithmState<Real> & algo_state)379 void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, 380 BoundConstraint<Real> &con, 381 AlgorithmState<Real> &algo_state ) {} 382 383 /** \brief Print iterate header. 384 */ printHeader(void) const385 std::string printHeader( void ) const { 386 std::stringstream hist; 387 hist << " "; 388 hist << std::setw(6) << std::left << "iter"; 389 hist << std::setw(15) << std::left << "fval"; 390 hist << std::setw(15) << std::left << "cnorm"; 391 hist << std::setw(15) << std::left << "gLnorm"; 392 hist << std::setw(15) << std::left << "snorm"; 393 hist << std::setw(10) << std::left << "delta"; 394 hist << std::setw(10) << std::left << "nnorm"; 395 hist << std::setw(10) << std::left << "tnorm"; 396 hist << std::setw(8) << std::left << "#fval"; 397 hist << std::setw(8) << std::left << "#grad"; 398 hist << std::setw(8) << std::left << "iterCG"; 399 hist << std::setw(8) << std::left << "flagCG"; 400 hist << std::setw(8) << std::left << "accept"; 401 hist << std::setw(8) << std::left << "linsys"; 402 hist << "\n"; 403 return hist.str(); 404 } 405 printName(void) const406 std::string printName( void ) const { 407 std::stringstream hist; 408 hist << "\n" << " Composite-step trust-region solver"; 409 hist << "\n"; 410 return hist.str(); 411 } 412 413 /** \brief Print iterate status. 414 */ print(AlgorithmState<Real> & algo_state,bool pHeader=false) const415 std::string print( AlgorithmState<Real> & algo_state, bool pHeader = false ) const { 416 //const ROL::Ptr<const StepState<Real> >& step_state = Step<Real>::getStepState(); 417 418 std::stringstream hist; 419 hist << std::scientific << std::setprecision(6); 420 if ( algo_state.iter == 0 ) { 421 hist << printName(); 422 } 423 if ( pHeader ) { 424 hist << printHeader(); 425 } 426 if ( algo_state.iter == 0 ) { 427 hist << " "; 428 hist << std::setw(6) << std::left << algo_state.iter; 429 hist << std::setw(15) << std::left << algo_state.value; 430 hist << std::setw(15) << std::left << algo_state.cnorm; 431 hist << std::setw(15) << std::left << algo_state.gnorm; 432 hist << "\n"; 433 } 434 else { 435 hist << " "; 436 hist << std::setw(6) << std::left << algo_state.iter; 437 hist << std::setw(15) << std::left << algo_state.value; 438 hist << std::setw(15) << std::left << algo_state.cnorm; 439 hist << std::setw(15) << std::left << algo_state.gnorm; 440 hist << std::setw(15) << std::left << algo_state.snorm; 441 hist << std::scientific << std::setprecision(2); 442 hist << std::setw(10) << std::left << Delta_; 443 hist << std::setw(10) << std::left << nnorm_; 444 hist << std::setw(10) << std::left << tnorm_; 445 hist << std::scientific << std::setprecision(6); 446 hist << std::setw(8) << std::left << algo_state.nfval; 447 hist << std::setw(8) << std::left << algo_state.ngrad; 448 hist << std::setw(8) << std::left << iterCG_; 449 hist << std::setw(8) << std::left << flagCG_; 450 hist << std::setw(8) << std::left << flagAC_; 451 hist << std::left << totalCallLS_ << "/" << totalIterLS_; 452 hist << "\n"; 453 } 454 return hist.str(); 455 } 456 457 /** \brief Compute Lagrange multipliers by solving the least-squares 458 problem minimizing the gradient of the Lagrangian, via the 459 augmented system formulation. 460 461 @param[out] l is the updated Lagrange multiplier; a dual constraint-space vector 462 @param[in] x is the current iterate; an optimization-space vector 463 @param[in] gf is the gradient of the objective function; a dual optimization-space vector 464 @param[in] con is the equality constraint object 465 466 On return ... fill the blanks. 467 */ computeLagrangeMultiplier(Vector<Real> & l,const Vector<Real> & x,const Vector<Real> & gf,Constraint<Real> & con)468 void computeLagrangeMultiplier(Vector<Real> &l, const Vector<Real> &x, const Vector<Real> &gf, Constraint<Real> &con) { 469 470 Real one(1); 471 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); 472 std::vector<Real> augiters; 473 474 if (infoLM_) { 475 std::stringstream hist; 476 hist << "\n Lagrange multiplier step\n"; 477 std::cout << hist.str(); 478 } 479 480 /* Apply adjoint of constraint Jacobian to current multiplier. */ 481 ROL::Ptr<Vector<Real> > ajl = gvec_->clone(); 482 con.applyAdjointJacobian(*ajl, l, x, zerotol); 483 484 /* Form right-hand side of the augmented system. */ 485 ROL::Ptr<Vector<Real> > b1 = gvec_->clone(); 486 ROL::Ptr<Vector<Real> > b2 = cvec_->clone(); 487 // b1 is the negative gradient of the Lagrangian 488 b1->set(gf); b1->plus(*ajl); b1->scale(-one); 489 // b2 is zero 490 b2->zero(); 491 492 /* Declare left-hand side of augmented system. */ 493 ROL::Ptr<Vector<Real> > v1 = xvec_->clone(); 494 ROL::Ptr<Vector<Real> > v2 = lvec_->clone(); 495 496 /* Compute linear solver tolerance. */ 497 Real b1norm = b1->norm(); 498 Real tol = setTolOSS(lmhtol_*b1norm); 499 500 /* Solve augmented system. */ 501 augiters = con.solveAugmentedSystem(*v1, *v2, *b1, *b2, x, tol); 502 totalCallLS_++; 503 totalIterLS_ = totalIterLS_ + augiters.size(); 504 printInfoLS(augiters); 505 506 /* Return updated Lagrange multiplier. */ 507 // v2 is the multiplier update 508 l.plus(*v2); 509 510 } // computeLagrangeMultiplier 511 512 513 /** \brief Compute quasi-normal step by minimizing the norm of 514 the linearized constraint. 515 516 Compute an approximate solution of the problem 517 \f[ 518 \begin{array}{rl} 519 \min_{n} & \|c'(x_k)n + c(x_k)\|^2_{\mathcal{X}} \\ 520 \mbox{subject to} & \|n\|_{\mathcal{X}} \le \delta 521 \end{array} 522 \f] 523 The approximate solution is computed using Powell's dogleg 524 method. The dogleg path is computed using the Cauchy point and 525 a full Newton step. The path's intersection with the trust-region 526 constraint gives the quasi-normal step. 527 528 @param[out] n is the quasi-normal step; an optimization-space vector 529 @param[in] c is the value of equality constraints; a constraint-space vector 530 @param[in] x is the current iterate; an optimization-space vector 531 @param[in] delta is the trust-region radius for the quasi-normal step 532 @param[in] con is the equality constraint object 533 534 */ computeQuasinormalStep(Vector<Real> & n,const Vector<Real> & c,const Vector<Real> & x,Real delta,Constraint<Real> & con)535 void computeQuasinormalStep(Vector<Real> &n, const Vector<Real> &c, const Vector<Real> &x, Real delta, Constraint<Real> &con) { 536 537 if (infoQN_) { 538 std::stringstream hist; 539 hist << "\n Quasi-normal step\n"; 540 std::cout << hist.str(); 541 } 542 543 Real zero(0); 544 Real one(1); 545 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); //zero; 546 std::vector<Real> augiters; 547 548 /* Compute Cauchy step nCP. */ 549 ROL::Ptr<Vector<Real> > nCP = xvec_->clone(); 550 ROL::Ptr<Vector<Real> > nCPdual = gvec_->clone(); 551 ROL::Ptr<Vector<Real> > nN = xvec_->clone(); 552 ROL::Ptr<Vector<Real> > ctemp = cvec_->clone(); 553 ROL::Ptr<Vector<Real> > dualc0 = lvec_->clone(); 554 dualc0->set(c.dual()); 555 con.applyAdjointJacobian(*nCPdual, *dualc0, x, zerotol); 556 nCP->set(nCPdual->dual()); 557 con.applyJacobian(*ctemp, *nCP, x, zerotol); 558 559 Real normsquare_ctemp = ctemp->dot(*ctemp); 560 if (normsquare_ctemp != zero) { 561 nCP->scale( -(nCP->dot(*nCP))/normsquare_ctemp ); 562 } 563 564 /* If the Cauchy step nCP is outside the trust region, 565 return the scaled Cauchy step. */ 566 Real norm_nCP = nCP->norm(); 567 if (norm_nCP >= delta) { 568 n.set(*nCP); 569 n.scale( delta/norm_nCP ); 570 if (infoQN_) { 571 std::stringstream hist; 572 hist << " taking partial Cauchy step\n"; 573 std::cout << hist.str(); 574 } 575 return; 576 } 577 578 /* Compute 'Newton' step, for example, by solving a problem 579 related to finding the minimum norm solution of min || c(x_k)*s + c ||^2. */ 580 // Compute tolerance for linear solver. 581 con.applyJacobian(*ctemp, *nCP, x, zerotol); 582 ctemp->plus(c); 583 Real tol = setTolOSS(qntol_*ctemp->norm()); 584 // Form right-hand side. 585 ctemp->scale(-one); 586 nCPdual->set(nCP->dual()); 587 nCPdual->scale(-one); 588 // Declare left-hand side of augmented system. 589 ROL::Ptr<Vector<Real> > dn = xvec_->clone(); 590 ROL::Ptr<Vector<Real> > y = lvec_->clone(); 591 // Solve augmented system. 592 augiters = con.solveAugmentedSystem(*dn, *y, *nCPdual, *ctemp, x, tol); 593 totalCallLS_++; 594 totalIterLS_ = totalIterLS_ + augiters.size(); 595 printInfoLS(augiters); 596 597 nN->set(*dn); 598 nN->plus(*nCP); 599 600 /* Either take full or partial Newton step, depending on 601 the trust-region constraint. */ 602 Real norm_nN = nN->norm(); 603 if (norm_nN <= delta) { 604 // Take full feasibility step. 605 n.set(*nN); 606 if (infoQN_) { 607 std::stringstream hist; 608 hist << " taking full Newton step\n"; 609 std::cout << hist.str(); 610 } 611 } 612 else { 613 // Take convex combination n = nCP+tau*(nN-nCP), 614 // so that ||n|| = delta. In other words, solve 615 // scalar quadratic equation: ||nCP+tau*(nN-nCP)||^2 = delta^2. 616 Real aa = dn->dot(*dn); 617 Real bb = dn->dot(*nCP); 618 Real cc = norm_nCP*norm_nCP - delta*delta; 619 Real tau = (-bb+sqrt(bb*bb-aa*cc))/aa; 620 n.set(*nCP); 621 n.axpy(tau, *dn); 622 if (infoQN_) { 623 std::stringstream hist; 624 hist << " taking dogleg step\n"; 625 std::cout << hist.str(); 626 } 627 } 628 629 } // computeQuasinormalStep 630 631 632 /** \brief Solve tangential subproblem. 633 634 @param[out] t is the solution of the tangential subproblem; an optimization-space vector 635 @param[out] tCP is the Cauchy point for the tangential subproblem; an optimization-space vector 636 @param[out] Wg is the dual of the projected gradient of the Lagrangian; an optimization-space vector 637 @param[in] x is the current iterate; an optimization-space vector 638 @param[in] g is the gradient of the Lagrangian; a dual optimization-space vector 639 @param[in] n is the quasi-normal step; an optimization-space vector 640 @param[in] l is the Lagrange multiplier; a dual constraint-space vector 641 @param[in] delta is the trust-region radius for the tangential subproblem 642 @param[in] con is the equality constraint object 643 644 */ solveTangentialSubproblem(Vector<Real> & t,Vector<Real> & tCP,Vector<Real> & Wg,const Vector<Real> & x,const Vector<Real> & g,const Vector<Real> & n,const Vector<Real> & l,Real delta,Objective<Real> & obj,Constraint<Real> & con)645 void solveTangentialSubproblem(Vector<Real> &t, Vector<Real> &tCP, Vector<Real> &Wg, 646 const Vector<Real> &x, const Vector<Real> &g, const Vector<Real> &n, const Vector<Real> &l, 647 Real delta, Objective<Real> &obj, Constraint<Real> &con) { 648 649 /* Initialization of the CG step. */ 650 bool orthocheck = true; // set to true if want to check orthogonality 651 // of Wr and r, otherwise set to false 652 Real tol_ortho = 0.5; // orthogonality measure; represets a bound on norm( \hat{S}, 2), where 653 // \hat{S} is defined in Heinkenschloss/Ridzal., "A Matrix-Free Trust-Region SQP Method" 654 Real S_max = 1.0; // another orthogonality measure; norm(S) needs to be bounded by 655 // a modest constant; norm(S) is small if the approximation of 656 // the null space projector is good 657 Real zero(0); 658 Real one(1); 659 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); 660 std::vector<Real> augiters; 661 iterCG_ = 1; 662 flagCG_ = 0; 663 t.zero(); 664 tCP.zero(); 665 ROL::Ptr<Vector<Real> > r = gvec_->clone(); 666 ROL::Ptr<Vector<Real> > pdesc = xvec_->clone(); 667 ROL::Ptr<Vector<Real> > tprev = xvec_->clone(); 668 ROL::Ptr<Vector<Real> > Wr = xvec_->clone(); 669 ROL::Ptr<Vector<Real> > Hp = gvec_->clone(); 670 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone(); 671 ROL::Ptr<Vector<Real> > gtemp = gvec_->clone(); 672 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone(); 673 ROL::Ptr<Vector<Real> > czero = cvec_->clone(); 674 czero->zero(); 675 r->set(g); 676 obj.hessVec(*gtemp, n, x, zerotol); 677 r->plus(*gtemp); 678 if (useConHess_) { 679 con.applyAdjointHessian(*gtemp, l, n, x, zerotol); 680 r->plus(*gtemp); 681 } 682 Real normg = r->norm(); 683 Real normWg = zero; 684 Real pHp = zero; 685 Real rp = zero; 686 Real alpha = zero; 687 Real normp = zero; 688 Real normr = zero; 689 Real normt = zero; 690 std::vector<Real> normWr(maxiterCG_+1, zero); 691 692 std::vector<ROL::Ptr<Vector<Real > > > p; // stores search directions 693 std::vector<ROL::Ptr<Vector<Real > > > Hps; // stores duals of hessvec's applied to p's 694 std::vector<ROL::Ptr<Vector<Real > > > rs; // stores duals of residuals 695 std::vector<ROL::Ptr<Vector<Real > > > Wrs; // stores duals of projected residuals 696 697 Real rptol(1e-12); 698 699 if (infoTS_) { 700 std::stringstream hist; 701 hist << "\n Tangential subproblem\n"; 702 hist << std::setw(6) << std::right << "iter" << std::setw(18) << "||Wr||/||Wr0||" << std::setw(15) << "||s||"; 703 hist << std::setw(15) << "delta" << std::setw(15) << "||c'(x)s||" << "\n"; 704 std::cout << hist.str(); 705 } 706 707 if (normg == 0) { 708 if (infoTS_) { 709 std::stringstream hist; 710 hist << " >>> Tangential subproblem: Initial gradient is zero! \n"; 711 std::cout << hist.str(); 712 } 713 iterCG_ = 0; Wg.zero(); flagCG_ = 0; 714 return; 715 } 716 717 /* Start CG loop. */ 718 while (iterCG_ < maxiterCG_) { 719 720 // Store tangential Cauchy point (which is the current iterate in the second iteration). 721 if (iterCG_ == 2) { 722 tCP.set(t); 723 } 724 725 // Compute (inexact) projection W*r. 726 if (iterCG_ == 1) { 727 // Solve augmented system. 728 Real tol = setTolOSS(pgtol_); 729 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol); 730 totalCallLS_++; 731 totalIterLS_ = totalIterLS_ + augiters.size(); 732 printInfoLS(augiters); 733 734 Wg.set(*Wr); 735 normWg = Wg.norm(); 736 if (orthocheck) { 737 Wrs.push_back(xvec_->clone()); 738 (Wrs[iterCG_-1])->set(*Wr); 739 } 740 // Check if done (small initial projected residual). 741 if (normWg == zero) { 742 flagCG_ = 0; 743 iterCG_--; 744 if (infoTS_) { 745 std::stringstream hist; 746 hist << " Initial projected residual is close to zero! \n"; 747 std::cout << hist.str(); 748 } 749 return; 750 } 751 // Set first residual to projected gradient. 752 // change r->set(Wg); 753 r->set(Wg.dual()); 754 if (orthocheck) { 755 rs.push_back(xvec_->clone()); 756 // change (rs[0])->set(*r); 757 (rs[0])->set(r->dual()); 758 } 759 } 760 else { 761 // Solve augmented system. 762 Real tol = setTolOSS(projtol_); 763 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol); 764 totalCallLS_++; 765 totalIterLS_ = totalIterLS_ + augiters.size(); 766 printInfoLS(augiters); 767 768 if (orthocheck) { 769 Wrs.push_back(xvec_->clone()); 770 (Wrs[iterCG_-1])->set(*Wr); 771 } 772 } 773 774 normWr[iterCG_-1] = Wr->norm(); 775 776 if (infoTS_) { 777 ROL::Ptr<Vector<Real> > ct = cvec_->clone(); 778 con.applyJacobian(*ct, t, x, zerotol); 779 Real linc = ct->norm(); 780 std::stringstream hist; 781 hist << std::scientific << std::setprecision(6); 782 hist << std::setw(6) << std::right << iterCG_-1 << std::setw(18) << normWr[iterCG_-1]/normWg << std::setw(15) << t.norm(); 783 hist << std::setw(15) << delta << std::setw(15) << linc << "\n"; 784 std::cout << hist.str(); 785 } 786 787 // Check if done (small relative residual). 788 if (normWr[iterCG_-1]/normWg < tolCG_) { 789 flagCG_ = 0; 790 iterCG_ = iterCG_-1; 791 if (infoTS_) { 792 std::stringstream hist; 793 hist << " || W(g + H*(n+s)) || <= cgtol*|| W(g + H*n)|| \n"; 794 std::cout << hist.str(); 795 } 796 return; 797 } 798 799 // Check nonorthogonality, one-norm of (WR*R/diag^2 - I) 800 if (orthocheck) { 801 ROL::LA::Matrix<Real> Wrr(iterCG_,iterCG_); // holds matrix Wrs'*rs 802 ROL::LA::Matrix<Real> T(iterCG_,iterCG_); // holds matrix T=(1/diag)*Wrs'*rs*(1/diag) 803 ROL::LA::Matrix<Real> Tm1(iterCG_,iterCG_); // holds matrix Tm1=T-I 804 for (int i=0; i<iterCG_; i++) { 805 for (int j=0; j<iterCG_; j++) { 806 Wrr(i,j) = (Wrs[i])->dot(*rs[j]); 807 T(i,j) = Wrr(i,j)/(normWr[i]*normWr[j]); 808 Tm1(i,j) = T(i,j); 809 if (i==j) { 810 Tm1(i,j) = Tm1(i,j) - one; 811 } 812 } 813 } 814 if (Tm1.normOne() >= tol_ortho) { 815 ROL::LAPACK<int,Real> lapack; 816 std::vector<int> ipiv(iterCG_); 817 int info; 818 std::vector<Real> work(3*iterCG_); 819 // compute inverse of T 820 lapack.GETRF(iterCG_, iterCG_, T.values(), T.stride(), &ipiv[0], &info); 821 lapack.GETRI(iterCG_, T.values(), T.stride(), &ipiv[0], &work[0], 3*iterCG_, &info); 822 Tm1 = T; 823 for (int i=0; i<iterCG_; i++) { 824 Tm1(i,i) = Tm1(i,i) - one; 825 } 826 if (Tm1.normOne() > S_max) { 827 flagCG_ = 4; 828 if (infoTS_) { 829 std::stringstream hist; 830 hist << " large nonorthogonality in W(R)'*R detected \n"; 831 std::cout << hist.str(); 832 } 833 return; 834 } 835 } 836 } 837 838 // Full orthogonalization. 839 p.push_back(xvec_->clone()); 840 (p[iterCG_-1])->set(*Wr); 841 (p[iterCG_-1])->scale(-one); 842 for (int j=1; j<iterCG_; j++) { 843 Real scal = (p[iterCG_-1])->dot(*(Hps[j-1])) / (p[j-1])->dot(*(Hps[j-1])); 844 ROL::Ptr<Vector<Real> > pj = xvec_->clone(); 845 pj->set(*p[j-1]); 846 pj->scale(-scal); 847 (p[iterCG_-1])->plus(*pj); 848 } 849 850 // change Hps.push_back(gvec_->clone()); 851 Hps.push_back(xvec_->clone()); 852 // change obj.hessVec(*(Hps[iterCG_-1]), *(p[iterCG_-1]), x, zerotol); 853 obj.hessVec(*Hp, *(p[iterCG_-1]), x, zerotol); 854 if (useConHess_) { 855 con.applyAdjointHessian(*gtemp, l, *(p[iterCG_-1]), x, zerotol); 856 // change (Hps[iterCG_-1])->plus(*gtemp); 857 Hp->plus(*gtemp); 858 } 859 // "Preconditioning" step. 860 (Hps[iterCG_-1])->set(Hp->dual()); 861 862 pHp = (p[iterCG_-1])->dot(*(Hps[iterCG_-1])); 863 // change rp = (p[iterCG_-1])->dot(*r); 864 rp = (p[iterCG_-1])->dot(*(rs[iterCG_-1])); 865 866 normp = (p[iterCG_-1])->norm(); 867 normr = r->norm(); 868 869 // Negative curvature stopping condition. 870 if (pHp <= 0) { 871 pdesc->set(*(p[iterCG_-1])); // p is the descent direction 872 if ((std::abs(rp) >= rptol*normp*normr) && (sgn(rp) == 1)) { 873 pdesc->scale(-one); // -p is the descent direction 874 } 875 flagCG_ = 2; 876 Real a = pdesc->dot(*pdesc); 877 Real b = pdesc->dot(t); 878 Real c = t.dot(t) - delta*delta; 879 // Positive root of a*theta^2 + 2*b*theta + c = 0. 880 Real theta = (-b + std::sqrt(b*b - a*c)) / a; 881 xtemp->set(*(p[iterCG_-1])); 882 xtemp->scale(theta); 883 t.plus(*xtemp); 884 // Store as tangential Cauchy point if terminating in first iteration. 885 if (iterCG_ == 1) { 886 tCP.set(t); 887 } 888 if (infoTS_) { 889 std::stringstream hist; 890 hist << " negative curvature detected \n"; 891 std::cout << hist.str(); 892 } 893 return; 894 } 895 896 // Want to enforce nonzero alpha's. 897 if (std::abs(rp) < rptol*normp*normr) { 898 flagCG_ = 5; 899 if (infoTS_) { 900 std::stringstream hist; 901 hist << " Zero alpha due to inexactness. \n"; 902 std::cout << hist.str(); 903 } 904 return; 905 } 906 907 alpha = - rp/pHp; 908 909 // Iterate update. 910 tprev->set(t); 911 xtemp->set(*(p[iterCG_-1])); 912 xtemp->scale(alpha); 913 t.plus(*xtemp); 914 915 // Trust-region stopping condition. 916 normt = t.norm(); 917 if (normt >= delta) { 918 pdesc->set(*(p[iterCG_-1])); // p is the descent direction 919 if (sgn(rp) == 1) { 920 pdesc->scale(-one); // -p is the descent direction 921 } 922 Real a = pdesc->dot(*pdesc); 923 Real b = pdesc->dot(*tprev); 924 Real c = tprev->dot(*tprev) - delta*delta; 925 // Positive root of a*theta^2 + 2*b*theta + c = 0. 926 Real theta = (-b + std::sqrt(b*b - a*c)) / a; 927 xtemp->set(*(p[iterCG_-1])); 928 xtemp->scale(theta); 929 t.set(*tprev); 930 t.plus(*xtemp); 931 // Store as tangential Cauchy point if terminating in first iteration. 932 if (iterCG_ == 1) { 933 tCP.set(t); 934 } 935 flagCG_ = 3; 936 if (infoTS_) { 937 std::stringstream hist; 938 hist << " trust-region condition active \n"; 939 std::cout << hist.str(); 940 } 941 return; 942 } 943 944 // Residual update. 945 xtemp->set(*(Hps[iterCG_-1])); 946 xtemp->scale(alpha); 947 // change r->plus(*gtemp); 948 r->plus(xtemp->dual()); 949 if (orthocheck) { 950 // change rs.push_back(gvec_->clone()); 951 rs.push_back(xvec_->clone()); 952 // change (rs[iterCG_])->set(*r); 953 (rs[iterCG_])->set(r->dual()); 954 } 955 956 iterCG_++; 957 958 } // while (iterCG_ < maxiterCG_) 959 960 flagCG_ = 1; 961 if (infoTS_) { 962 std::stringstream hist; 963 hist << " maximum number of iterations reached \n"; 964 std::cout << hist.str(); 965 } 966 967 } // solveTangentialSubproblem 968 969 970 /** \brief Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence. 971 */ accept(Vector<Real> & s,Vector<Real> & n,Vector<Real> & t,Real f_new,Vector<Real> & c_new,Vector<Real> & gf_new,Vector<Real> & l_new,Vector<Real> & g_new,const Vector<Real> & x,const Vector<Real> & l,Real f,const Vector<Real> & gf,const Vector<Real> & c,const Vector<Real> & g,Vector<Real> & tCP,Vector<Real> & Wg,Objective<Real> & obj,Constraint<Real> & con,AlgorithmState<Real> & algo_state)972 void accept(Vector<Real> &s, Vector<Real> &n, Vector<Real> &t, Real f_new, Vector<Real> &c_new, 973 Vector<Real> &gf_new, Vector<Real> &l_new, Vector<Real> &g_new, 974 const Vector<Real> &x, const Vector<Real> &l, Real f, const Vector<Real> &gf, const Vector<Real> &c, 975 const Vector<Real> &g, Vector<Real> &tCP, Vector<Real> &Wg, 976 Objective<Real> &obj, Constraint<Real> &con, AlgorithmState<Real> &algo_state) { 977 978 Real beta = 1e-8; // predicted reduction parameter 979 Real tol_red_tang = 1e-3; // internal reduction factor for tangtol 980 Real tol_red_all = 1e-1; // internal reduction factor for qntol, lmhtol, pgtol, projtol, tangtol 981 //bool glob_refine = true; // true - if subsolver tolerances are adjusted in this routine, keep adjusted values globally 982 // false - if subsolver tolerances are adjusted in this routine, discard adjusted values 983 Real tol_fdiff = 1e-12; // relative objective function difference for ared computation 984 int ct_max = 10; // maximum number of globalization tries 985 Real mintol = 1e-16; // smallest projection tolerance value 986 987 // Determines max value of |rpred|/pred. 988 Real rpred_over_pred = 0.5*(1-eta_); 989 990 if (infoAC_) { 991 std::stringstream hist; 992 hist << "\n Composite step acceptance\n"; 993 std::cout << hist.str(); 994 } 995 996 Real zero = 0.0; 997 Real one = 1.0; 998 Real two = 2.0; 999 Real half = one/two; 1000 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); 1001 std::vector<Real> augiters; 1002 1003 Real pred = zero; 1004 Real ared = zero; 1005 Real rpred = zero; 1006 Real part_pred = zero; 1007 Real linc_preproj = zero; 1008 Real linc_postproj = zero; 1009 Real tangtol_start = zero; 1010 Real tangtol = tangtol_; 1011 //Real projtol = projtol_; 1012 bool flag = false; 1013 int num_proj = 0; 1014 bool try_tCP = false; 1015 Real fdiff = zero; 1016 1017 ROL::Ptr<Vector<Real> > xtrial = xvec_->clone(); 1018 ROL::Ptr<Vector<Real> > Jl = gvec_->clone(); 1019 ROL::Ptr<Vector<Real> > gfJl = gvec_->clone(); 1020 ROL::Ptr<Vector<Real> > Jnc = cvec_->clone(); 1021 ROL::Ptr<Vector<Real> > t_orig = xvec_->clone(); 1022 ROL::Ptr<Vector<Real> > t_dual = gvec_->clone(); 1023 ROL::Ptr<Vector<Real> > Jt_orig = cvec_->clone(); 1024 ROL::Ptr<Vector<Real> > t_m_tCP = xvec_->clone(); 1025 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone(); 1026 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone(); 1027 ROL::Ptr<Vector<Real> > rt = cvec_->clone(); 1028 ROL::Ptr<Vector<Real> > Hn = gvec_->clone(); 1029 ROL::Ptr<Vector<Real> > Hto = gvec_->clone(); 1030 ROL::Ptr<Vector<Real> > cxxvec = gvec_->clone(); 1031 ROL::Ptr<Vector<Real> > czero = cvec_->clone(); 1032 czero->zero(); 1033 Real Jnc_normsquared = zero; 1034 Real c_normsquared = zero; 1035 1036 // Compute and store some quantities for later use. Necessary 1037 // because of the function and constraint updates below. 1038 con.applyAdjointJacobian(*Jl, l, x, zerotol); 1039 con.applyJacobian(*Jnc, n, x, zerotol); 1040 Jnc->plus(c); 1041 Jnc_normsquared = Jnc->dot(*Jnc); 1042 c_normsquared = c.dot(c); 1043 1044 for (int ct=0; ct<ct_max; ct++) { 1045 1046 try_tCP = true; 1047 t_m_tCP->set(t); 1048 t_m_tCP->scale(-one); 1049 t_m_tCP->plus(tCP); 1050 if (t_m_tCP->norm() == zero) { 1051 try_tCP = false; 1052 } 1053 1054 t_orig->set(t); 1055 con.applyJacobian(*Jt_orig, *t_orig, x, zerotol); 1056 linc_preproj = Jt_orig->norm(); 1057 pred = one; 1058 rpred = two*rpred_over_pred*pred; 1059 flag = false; 1060 num_proj = 1; 1061 tangtol_start = tangtol; 1062 1063 while (std::abs(rpred)/pred > rpred_over_pred) { 1064 // Compute projected tangential step. 1065 if (flag) { 1066 tangtol = tol_red_tang*tangtol; 1067 num_proj++; 1068 if (tangtol < mintol) { 1069 if (infoAC_) { 1070 std::stringstream hist; 1071 hist << "\n The projection of the tangential step cannot be done with sufficient precision.\n"; 1072 hist << " Is the quasi-normal step very small? Continuing with no global convergence guarantees.\n"; 1073 std::cout << hist.str(); 1074 } 1075 break; 1076 } 1077 } 1078 // Solve augmented system. 1079 Real tol = setTolOSS(tangtol); 1080 // change augiters = con.solveAugmentedSystem(t, *ltemp, *t_orig, *czero, x, tol); 1081 t_dual->set(t_orig->dual()); 1082 augiters = con.solveAugmentedSystem(t, *ltemp, *t_dual, *czero, x, tol); 1083 totalCallLS_++; 1084 totalIterLS_ = totalIterLS_ + augiters.size(); 1085 printInfoLS(augiters); 1086 totalProj_++; 1087 con.applyJacobian(*rt, t, x, zerotol); 1088 linc_postproj = rt->norm(); 1089 1090 // Compute composite step. 1091 s.set(t); 1092 s.plus(n); 1093 1094 // Compute some quantities before updating the objective and the constraint. 1095 obj.hessVec(*Hn, n, x, zerotol); 1096 if (useConHess_) { 1097 con.applyAdjointHessian(*cxxvec, l, n, x, zerotol); 1098 Hn->plus(*cxxvec); 1099 } 1100 obj.hessVec(*Hto, *t_orig, x, zerotol); 1101 if (useConHess_) { 1102 con.applyAdjointHessian(*cxxvec, l, *t_orig, x, zerotol); 1103 Hto->plus(*cxxvec); 1104 } 1105 1106 // Compute objective, constraint, etc. values at the trial point. 1107 xtrial->set(x); 1108 xtrial->plus(s); 1109 obj.update(*xtrial,false,algo_state.iter); 1110 con.update(*xtrial,false,algo_state.iter); 1111 f_new = obj.value(*xtrial, zerotol); 1112 obj.gradient(gf_new, *xtrial, zerotol); 1113 con.value(c_new, *xtrial, zerotol); 1114 l_new.set(l); 1115 computeLagrangeMultiplier(l_new, *xtrial, gf_new, con); 1116 1117 // Penalty parameter update. 1118 part_pred = - Wg.dot(*t_orig); 1119 gfJl->set(gf); 1120 gfJl->plus(*Jl); 1121 // change part_pred -= gfJl->dot(n); 1122 part_pred -= n.dot(gfJl->dual()); 1123 // change part_pred -= half*Hn->dot(n); 1124 part_pred -= half*n.dot(Hn->dual()); 1125 // change part_pred -= half*Hto->dot(*t_orig); 1126 part_pred -= half*t_orig->dot(Hto->dual()); 1127 ltemp->set(l_new); 1128 ltemp->axpy(-one, l); 1129 // change part_pred -= Jnc->dot(*ltemp); 1130 part_pred -= Jnc->dot(ltemp->dual()); 1131 1132 if ( part_pred < -half*penalty_*(c_normsquared-Jnc_normsquared) ) { 1133 penalty_ = ( -two * part_pred / (c_normsquared-Jnc_normsquared) ) + beta; 1134 } 1135 1136 pred = part_pred + penalty_*(c_normsquared-Jnc_normsquared); 1137 1138 // Computation of rpred. 1139 // change rpred = - ltemp->dot(*rt) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc); 1140 rpred = - rt->dot(ltemp->dual()) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc); 1141 // change ROL::Ptr<Vector<Real> > lrt = lvec_->clone(); 1142 //lrt->set(*rt); 1143 //rpred = - ltemp->dot(*rt) - penalty_ * std::pow(rt->norm(), 2) - two * penalty_ * lrt->dot(*Jnc); 1144 flag = 1; 1145 1146 } // while (std::abs(rpred)/pred > rpred_over_pred) 1147 1148 tangtol = tangtol_start; 1149 1150 // Check if the solution of the tangential subproblem is 1151 // disproportionally large compared to total trial step. 1152 xtemp->set(n); 1153 xtemp->plus(t); 1154 if ( t_orig->norm()/xtemp->norm() < tntmax_ ) { 1155 break; 1156 } 1157 else { 1158 t_m_tCP->set(*t_orig); 1159 t_m_tCP->scale(-one); 1160 t_m_tCP->plus(tCP); 1161 if ((t_m_tCP->norm() > 0) && try_tCP) { 1162 if (infoAC_) { 1163 std::stringstream hist; 1164 hist << " ---> now trying tangential Cauchy point\n"; 1165 std::cout << hist.str(); 1166 } 1167 t.set(tCP); 1168 } 1169 else { 1170 if (infoAC_) { 1171 std::stringstream hist; 1172 hist << " ---> recomputing quasi-normal step and re-solving tangential subproblem\n"; 1173 std::cout << hist.str(); 1174 } 1175 totalRef_++; 1176 // Reset global quantities. 1177 obj.update(x); 1178 con.update(x); 1179 /*lmhtol = tol_red_all*lmhtol; 1180 qntol = tol_red_all*qntol; 1181 pgtol = tol_red_all*pgtol; 1182 projtol = tol_red_all*projtol; 1183 tangtol = tol_red_all*tangtol; 1184 if (glob_refine) { 1185 lmhtol_ = lmhtol; 1186 qntol_ = qntol; 1187 pgtol_ = pgtol; 1188 projtol_ = projtol; 1189 tangtol_ = tangtol; 1190 }*/ 1191 if (!tolOSSfixed_) { 1192 lmhtol_ *= tol_red_all; 1193 qntol_ *= tol_red_all; 1194 pgtol_ *= tol_red_all; 1195 projtol_ *= tol_red_all; 1196 tangtol_ *= tol_red_all; 1197 } 1198 // Recompute the quasi-normal step. 1199 computeQuasinormalStep(n, c, x, zeta_*Delta_, con); 1200 // Solve tangential subproblem. 1201 solveTangentialSubproblem(t, tCP, Wg, x, g, n, l, Delta_, obj, con); 1202 totalIterCG_ += iterCG_; 1203 if (flagCG_ == 1) { 1204 totalNegCurv_++; 1205 } 1206 } 1207 } // else w.r.t. ( t_orig->norm()/xtemp->norm() < tntmax ) 1208 1209 } // for (int ct=0; ct<ct_max; ct++) 1210 1211 // Compute actual reduction; 1212 fdiff = f - f_new; 1213 // Heuristic 1: If fdiff is very small compared to f, set it to 0, 1214 // in order to prevent machine precision issues. 1215 Real em24(1e-24); 1216 Real em14(1e-14); 1217 if (std::abs(fdiff / (f+em24)) < tol_fdiff) { 1218 fdiff = em14; 1219 } 1220 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(c.dot(c) - c_new.dot(c_new)); 1221 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(std::pow(c.norm(),2) - std::pow(c_new.norm(),2)); 1222 ared = fdiff + (c.dot(l.dual()) - c_new.dot(l_new.dual())) + penalty_*(c.dot(c) - c_new.dot(c_new)); 1223 1224 // Store actual and predicted reduction. 1225 ared_ = ared; 1226 pred_ = pred; 1227 1228 // Store step and vector norms. 1229 snorm_ = s.norm(); 1230 nnorm_ = n.norm(); 1231 tnorm_ = t.norm(); 1232 1233 // Print diagnostics. 1234 if (infoAC_) { 1235 std::stringstream hist; 1236 hist << "\n Trial step info ...\n"; 1237 hist << " n_norm = " << nnorm_ << "\n"; 1238 hist << " t_norm = " << tnorm_ << "\n"; 1239 hist << " s_norm = " << snorm_ << "\n"; 1240 hist << " xtrial_norm = " << xtrial->norm() << "\n"; 1241 hist << " f_old = " << f << "\n"; 1242 hist << " f_trial = " << f_new << "\n"; 1243 hist << " f_old-f_trial = " << f-f_new << "\n"; 1244 hist << " ||c_old|| = " << c.norm() << "\n"; 1245 hist << " ||c_trial|| = " << c_new.norm() << "\n"; 1246 hist << " ||Jac*t_preproj|| = " << linc_preproj << "\n"; 1247 hist << " ||Jac*t_postproj|| = " << linc_postproj << "\n"; 1248 hist << " ||t_tilde||/||t|| = " << t_orig->norm() / t.norm() << "\n"; 1249 hist << " ||t_tilde||/||n+t|| = " << t_orig->norm() / snorm_ << "\n"; 1250 hist << " # projections = " << num_proj << "\n"; 1251 hist << " penalty param = " << penalty_ << "\n"; 1252 hist << " ared = " << ared_ << "\n"; 1253 hist << " pred = " << pred_ << "\n"; 1254 hist << " ared/pred = " << ared_/pred_ << "\n"; 1255 std::cout << hist.str(); 1256 } 1257 1258 } // accept 1259 1260 }; // class CompositeStep 1261 1262 } // namespace ROL 1263 1264 #endif 1265