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_ScalarMinimizationLineSearch_H 45 #define ROL_ScalarMinimizationLineSearch_H 46 47 /** \class ROL::ScalarMinimizationLineSearch 48 \brief Implements line search methods that attempt to minimize the 49 scalar function \f$\phi(t) := f(x+ts)\f$. 50 */ 51 52 #include "ROL_LineSearch.hpp" 53 #include "ROL_BrentsScalarMinimization.hpp" 54 #include "ROL_BisectionScalarMinimization.hpp" 55 #include "ROL_GoldenSectionScalarMinimization.hpp" 56 #include "ROL_ScalarFunction.hpp" 57 #include "ROL_Bracketing.hpp" 58 59 namespace ROL { 60 61 template<class Real> 62 class ScalarMinimizationLineSearch : public LineSearch<Real> { 63 private: 64 ROL::Ptr<Vector<Real> > xnew_; 65 ROL::Ptr<Vector<Real> > g_; 66 ROL::Ptr<ScalarMinimization<Real> > sm_; 67 ROL::Ptr<Bracketing<Real> > br_; 68 ROL::Ptr<ScalarFunction<Real> > sf_; 69 70 ECurvatureCondition econd_; 71 Real c1_; 72 Real c2_; 73 Real c3_; 74 int max_nfval_; 75 76 class Phi : public ScalarFunction<Real> { 77 private: 78 const ROL::Ptr<Vector<Real> > xnew_; 79 const ROL::Ptr<Vector<Real> > g_; 80 const ROL::Ptr<const Vector<Real> > x_; 81 const ROL::Ptr<const Vector<Real> > s_; 82 const ROL::Ptr<Objective<Real> > obj_; 83 const ROL::Ptr<BoundConstraint<Real> > con_; 84 Real ftol_; updateIterate(Real alpha)85 void updateIterate(Real alpha) { 86 xnew_->set(*x_); 87 xnew_->axpy(alpha,*s_); 88 if ( con_->isActivated() ) { 89 con_->project(*xnew_); 90 } 91 } 92 public: Phi(const ROL::Ptr<Vector<Real>> & xnew,const ROL::Ptr<Vector<Real>> & g,const ROL::Ptr<const Vector<Real>> & x,const ROL::Ptr<const Vector<Real>> & s,const ROL::Ptr<Objective<Real>> & obj,const ROL::Ptr<BoundConstraint<Real>> & con)93 Phi(const ROL::Ptr<Vector<Real> > &xnew, 94 const ROL::Ptr<Vector<Real> > &g, 95 const ROL::Ptr<const Vector<Real> > &x, 96 const ROL::Ptr<const Vector<Real> > &s, 97 const ROL::Ptr<Objective<Real> > &obj, 98 const ROL::Ptr<BoundConstraint<Real> > &con) 99 : xnew_(xnew), g_(g), x_(x), s_(s), obj_(obj), con_(con), 100 ftol_(std::sqrt(ROL_EPSILON<Real>())) {} value(const Real alpha)101 Real value(const Real alpha) { 102 updateIterate(alpha); 103 obj_->update(*xnew_); 104 return obj_->value(*xnew_,ftol_); 105 } deriv(const Real alpha)106 Real deriv(const Real alpha) { 107 updateIterate(alpha); 108 obj_->update(*xnew_); 109 obj_->gradient(*g_,*xnew_,ftol_); 110 return s_->dot(g_->dual()); 111 } 112 }; 113 114 class LineSearchStatusTest : public ScalarMinimizationStatusTest<Real> { 115 private: 116 ROL::Ptr<ScalarFunction<Real> > phi_; 117 118 const Real f0_; 119 const Real g0_; 120 121 const Real c1_; 122 const Real c2_; 123 const Real c3_; 124 const int max_nfval_; 125 const ECurvatureCondition econd_; 126 127 128 public: LineSearchStatusTest(const Real f0,const Real g0,const Real c1,const Real c2,const Real c3,const int max_nfval,ECurvatureCondition econd,const ROL::Ptr<ScalarFunction<Real>> & phi)129 LineSearchStatusTest(const Real f0, const Real g0, 130 const Real c1, const Real c2, const Real c3, 131 const int max_nfval, ECurvatureCondition econd, 132 const ROL::Ptr<ScalarFunction<Real> > &phi) 133 : phi_(phi), f0_(f0), g0_(g0), c1_(c1), c2_(c2), c3_(c3), 134 max_nfval_(max_nfval), econd_(econd) {} 135 check(Real & x,Real & fx,Real & gx,int & nfval,int & ngval,const bool deriv=false)136 bool check(Real &x, Real &fx, Real &gx, 137 int &nfval, int &ngval, const bool deriv = false) { 138 Real one(1), two(2); 139 bool armijo = (fx <= f0_ + c1_*x*g0_); 140 // bool itcond = (nfval >= max_nfval_); 141 bool curvcond = false; 142 // if (armijo && !itcond) { 143 if (armijo) { 144 if (econd_ == CURVATURECONDITION_GOLDSTEIN) { 145 curvcond = (fx >= f0_ + (one-c1_)*x*g0_); 146 } 147 else if (econd_ == CURVATURECONDITION_NULL) { 148 curvcond = true; 149 } 150 else { 151 if (!deriv) { 152 gx = phi_->deriv(x); ngval++; 153 } 154 if (econd_ == CURVATURECONDITION_WOLFE) { 155 curvcond = (gx >= c2_*g0_); 156 } 157 else if (econd_ == CURVATURECONDITION_STRONGWOLFE) { 158 curvcond = (std::abs(gx) <= c2_*std::abs(g0_)); 159 } 160 else if (econd_ == CURVATURECONDITION_GENERALIZEDWOLFE) { 161 curvcond = (c2_*g0_ <= gx && gx <= -c3_*g0_); 162 } 163 else if (econd_ == CURVATURECONDITION_APPROXIMATEWOLFE) { 164 curvcond = (c2_*g0_ <= gx && gx <= (two*c1_ - one)*g0_); 165 } 166 } 167 } 168 //return (armijo && curvcond) || itcond; 169 return (armijo && curvcond); 170 } 171 }; 172 173 public: 174 // Constructor ScalarMinimizationLineSearch(ROL::ParameterList & parlist,const ROL::Ptr<ScalarMinimization<Real>> & sm=ROL::nullPtr,const ROL::Ptr<Bracketing<Real>> & br=ROL::nullPtr,const ROL::Ptr<ScalarFunction<Real>> & sf=ROL::nullPtr)175 ScalarMinimizationLineSearch( ROL::ParameterList &parlist, 176 const ROL::Ptr<ScalarMinimization<Real> > &sm = ROL::nullPtr, 177 const ROL::Ptr<Bracketing<Real> > &br = ROL::nullPtr, 178 const ROL::Ptr<ScalarFunction<Real> > &sf = ROL::nullPtr ) 179 : LineSearch<Real>(parlist) { 180 Real zero(0), p4(0.4), p6(0.6), p9(0.9), oem4(1.e-4), oem10(1.e-10), one(1); 181 ROL::ParameterList &list0 = parlist.sublist("Step").sublist("Line Search"); 182 ROL::ParameterList &list = list0.sublist("Line-Search Method"); 183 // Get Bracketing Method 184 if( br == ROL::nullPtr ) { 185 br_ = ROL::makePtr<Bracketing<Real>>(); 186 } 187 else { 188 br_ = br; 189 } 190 // Get ScalarMinimization Method 191 std::string type = list.get("Type","Brent's"); 192 Real tol = list.sublist(type).get("Tolerance",oem10); 193 int niter = list.sublist(type).get("Iteration Limit",1000); 194 ROL::ParameterList plist; 195 plist.sublist("Scalar Minimization").set("Type",type); 196 plist.sublist("Scalar Minimization").sublist(type).set("Tolerance",tol); 197 plist.sublist("Scalar Minimization").sublist(type).set("Iteration Limit",niter); 198 199 if( sm == ROL::nullPtr ) { // No user-provided ScalarMinimization object 200 201 if ( type == "Brent's" ) { 202 sm_ = ROL::makePtr<BrentsScalarMinimization<Real>>(plist); 203 } 204 else if ( type == "Bisection" ) { 205 sm_ = ROL::makePtr<BisectionScalarMinimization<Real>>(plist); 206 } 207 else if ( type == "Golden Section" ) { 208 sm_ = ROL::makePtr<GoldenSectionScalarMinimization<Real>>(plist); 209 } 210 else { 211 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument, 212 ">>> (ROL::ScalarMinimizationLineSearch): Undefined ScalarMinimization type!"); 213 } 214 } 215 else { 216 sm_ = sm; 217 } 218 219 sf_ = sf; 220 221 222 // Status test for line search 223 econd_ = StringToECurvatureCondition(list0.sublist("Curvature Condition").get("Type","Strong Wolfe Conditions")); 224 max_nfval_ = list0.get("Function Evaluation Limit",20); 225 c1_ = list0.get("Sufficient Decrease Tolerance",oem4); 226 c2_ = list0.sublist("Curvature Condition").get("General Parameter",p9); 227 c3_ = list0.sublist("Curvature Condition").get("Generalized Wolfe Parameter",p6); 228 // Check status test inputs 229 c1_ = ((c1_ < zero) ? oem4 : c1_); 230 c2_ = ((c2_ < zero) ? p9 : c2_); 231 c3_ = ((c3_ < zero) ? p9 : c3_); 232 if ( c2_ <= c1_ ) { 233 c1_ = oem4; 234 c2_ = p9; 235 } 236 EDescent edesc = StringToEDescent(list0.sublist("Descent Method").get("Type","Quasi-Newton Method")); 237 if ( edesc == DESCENT_NONLINEARCG ) { 238 c2_ = p4; 239 c3_ = std::min(one-c2_,c3_); 240 } 241 } 242 initialize(const Vector<Real> & x,const Vector<Real> & s,const Vector<Real> & g,Objective<Real> & obj,BoundConstraint<Real> & con)243 void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g, 244 Objective<Real> &obj, BoundConstraint<Real> &con ) { 245 LineSearch<Real>::initialize(x,s,g,obj,con); 246 xnew_ = x.clone(); 247 g_ = g.clone(); 248 } 249 250 // Find the minimum of phi(alpha) = f(x + alpha*s) using Brent's method run(Real & alpha,Real & fval,int & ls_neval,int & ls_ngrad,const Real & gs,const Vector<Real> & s,const Vector<Real> & x,Objective<Real> & obj,BoundConstraint<Real> & con)251 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad, 252 const Real &gs, const Vector<Real> &s, const Vector<Real> &x, 253 Objective<Real> &obj, BoundConstraint<Real> &con ) { 254 ls_neval = 0; ls_ngrad = 0; 255 256 // Get initial line search parameter 257 alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con); 258 259 // Build ScalarFunction and ScalarMinimizationStatusTest 260 ROL::Ptr<const Vector<Real> > x_ptr = ROL::makePtrFromRef(x); 261 ROL::Ptr<const Vector<Real> > s_ptr = ROL::makePtrFromRef(s); 262 ROL::Ptr<Objective<Real> > obj_ptr = ROL::makePtrFromRef(obj); 263 ROL::Ptr<BoundConstraint<Real> > bnd_ptr = ROL::makePtrFromRef(con); 264 265 266 ROL::Ptr<ScalarFunction<Real> > phi; 267 268 if( sf_ == ROL::nullPtr ) { 269 phi = ROL::makePtr<Phi>(xnew_,g_,x_ptr,s_ptr,obj_ptr,bnd_ptr); 270 } 271 else { 272 phi = sf_; 273 } 274 275 ROL::Ptr<ScalarMinimizationStatusTest<Real> > test 276 = ROL::makePtr<LineSearchStatusTest>(fval,gs,c1_,c2_,c3_,max_nfval_,econd_,phi); 277 278 // Run Bracketing 279 int nfval = 0, ngrad = 0; 280 Real A(0), fA = fval; 281 Real B = alpha, fB = phi->value(B); 282 br_->run(alpha,fval,A,fA,B,fB,nfval,ngrad,*phi,*test); 283 B = alpha; 284 ls_neval += nfval; ls_ngrad += ngrad; 285 286 // Run ScalarMinimization 287 nfval = 0, ngrad = 0; 288 sm_->run(fval, alpha, nfval, ngrad, *phi, A, B, *test); 289 ls_neval += nfval; ls_ngrad += ngrad; 290 291 LineSearch<Real>::setNextInitialAlpha(alpha); 292 } 293 }; 294 295 } 296 297 #endif 298