ROL
ROL_LineSearch.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_LINESEARCH_H
45 #define ROL_LINESEARCH_H
46 
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_ParameterList.hpp"
53 #include "ROL_Types.hpp"
54 #include "ROL_Vector.hpp"
55 #include "ROL_Objective.hpp"
56 #include "ROL_BoundConstraint.hpp"
57 #include "ROL_ScalarFunction.hpp"
58 
59 namespace ROL {
60 
61 template<class Real>
62 class LineSearch {
63 private:
64 
67 
68  bool useralpha_;
69  bool usePrevAlpha_; // Use the previous step's accepted alpha as an initial guess
70  Real alpha0_;
71  int maxit_;
72  Real c1_;
73  Real c2_;
74  Real c3_;
75  Real eps_;
76  Real fmin_; // smallest fval encountered
77  Real alphaMin_; // Alpha that yields the smallest fval encountered
78  bool acceptMin_; // Use smallest fval if sufficient decrease not satisfied
79  bool itcond_; // true if maximum function evaluations reached
80 
81  Teuchos::RCP<Vector<Real> > xtst_;
82  Teuchos::RCP<Vector<Real> > d_;
83  Teuchos::RCP<Vector<Real> > g_;
84  Teuchos::RCP<Vector<Real> > grad_;
85 // Teuchos::RCP<const Vector<Real> > grad_;
86 
87 public:
88 
89 
90  virtual ~LineSearch() {}
91 
92  // Constructor
93  LineSearch( Teuchos::ParameterList &parlist ) : eps_(0) {
94  Real one(1), p9(0.9), p6(0.6), p4(0.4), oem4(1.e-4), zero(0);
95  // Enumerations
96  edesc_ = StringToEDescent(parlist.sublist("Step").sublist("Line Search").sublist("Descent Method").get("Type","Quasi-Newton Method"));
97  econd_ = StringToECurvatureCondition(parlist.sublist("Step").sublist("Line Search").sublist("Curvature Condition").get("Type","Strong Wolfe Conditions"));
98  // Linesearch Parameters
99  alpha0_ = parlist.sublist("Step").sublist("Line Search").get("Initial Step Size",one);
100  useralpha_ = parlist.sublist("Step").sublist("Line Search").get("User Defined Initial Step Size",false);
101  usePrevAlpha_ = parlist.sublist("Step").sublist("Line Search").get("Use Previous Step Length as Initial Guess",false);
102  acceptMin_ = parlist.sublist("Step").sublist("Line Search").get("Accept Linesearch Minimizer",false);
103  maxit_ = parlist.sublist("Step").sublist("Line Search").get("Function Evaluation Limit",20);
104  c1_ = parlist.sublist("Step").sublist("Line Search").get("Sufficient Decrease Tolerance",oem4);
105  c2_ = parlist.sublist("Step").sublist("Line Search").sublist("Curvature Condition").get("General Parameter",p9);
106  c3_ = parlist.sublist("Step").sublist("Line Search").sublist("Curvature Condition").get("Generalized Wolfe Parameter",p6);
107 
108  fmin_ = std::numeric_limits<Real>::max();
109  alphaMin_ = 0;
110  itcond_ = false;
111 
112  c1_ = ((c1_ < zero) ? oem4 : c1_);
113  c2_ = ((c2_ < zero) ? p9 : c2_);
114  c3_ = ((c3_ < zero) ? p9 : c3_);
115  if ( c2_ <= c1_ ) {
116  c1_ = oem4;
117  c2_ = p9;
118  }
119  if ( edesc_ == DESCENT_NONLINEARCG ) {
120  c2_ = p4;
121  c3_ = std::min(one-c2_,c3_);
122  }
123  }
124 
125 
126  virtual void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
128  grad_ = g.clone();
129  xtst_ = x.clone();
130  d_ = s.clone();
131  g_ = g.clone();
132  }
133 
134  virtual void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
135  const Real &gs, const Vector<Real> &s, const Vector<Real> &x,
136  Objective<Real> &obj, BoundConstraint<Real> &con ) = 0;
137 
138  // Set epsilon for epsilon active sets
139  void setData(Real &eps, const Vector<Real> &g) {
140  eps_ = eps;
141  grad_->set(g);
142  }
143 
144  // use this function to modify alpha and fval if the maximum number of iterations
145  // are reached
146  void setMaxitUpdate(Real &alpha, Real &fnew, const Real &fold) {
147  // Use local minimizer
148  if( itcond_ && acceptMin_ ) {
149  alpha = alphaMin_;
150  fnew = fmin_;
151  }
152  // Take no step
153  else if(itcond_ && !acceptMin_) {
154  alpha = 0;
155  fnew = fold;
156  }
157  setNextInitialAlpha(alpha);
158  }
159 
160 
161 protected:
162  virtual bool status( const ELineSearch type, int &ls_neval, int &ls_ngrad, const Real alpha,
163  const Real fold, const Real sgold, const Real fnew,
164  const Vector<Real> &x, const Vector<Real> &s,
165  Objective<Real> &obj, BoundConstraint<Real> &con ) {
166  Real tol = std::sqrt(ROL_EPSILON<Real>()), one(1), two(2);
167 
168  // Check Armijo Condition
169  bool armijo = false;
170  if ( con.isActivated() ) {
171  Real gs(0);
172  if ( edesc_ == DESCENT_STEEPEST ) {
173  updateIterate(*d_,x,s,alpha,con);
174  d_->scale(-one);
175  d_->plus(x);
176  gs = -s.dot(*d_);
177  }
178  else {
179  d_->set(s);
180  d_->scale(-one);
181  con.pruneActive(*d_,grad_->dual(),x,eps_);
182  gs = alpha*(grad_)->dot(d_->dual());
183  d_->zero();
184  updateIterate(*d_,x,s,alpha,con);
185  d_->scale(-one);
186  d_->plus(x);
187  con.pruneInactive(*d_,grad_->dual(),x,eps_);
188  gs += d_->dot(grad_->dual());
189  }
190  if ( fnew <= fold - c1_*gs ) {
191  armijo = true;
192  }
193  }
194  else {
195  if ( fnew <= fold + c1_*alpha*sgold ) {
196  armijo = true;
197  }
198  }
199 
200  // Check Maximum Iteration
201  itcond_ = false;
202  if ( ls_neval >= maxit_ ) {
203  itcond_ = true;
204  }
205 
206  // Check Curvature Condition
207  bool curvcond = false;
208  if ( armijo && ((type != LINESEARCH_BACKTRACKING && type != LINESEARCH_CUBICINTERP) ||
209  (edesc_ == DESCENT_NONLINEARCG)) ) {
211  if (fnew >= fold + (one-c1_)*alpha*sgold) {
212  curvcond = true;
213  }
214  }
215  else if (econd_ == CURVATURECONDITION_NULL) {
216  curvcond = true;
217  }
218  else {
219  updateIterate(*xtst_,x,s,alpha,con);
220  obj.update(*xtst_);
221  obj.gradient(*g_,*xtst_,tol);
222  Real sgnew(0);
223  if ( con.isActivated() ) {
224  d_->set(s);
225  d_->scale(-alpha);
226  con.pruneActive(*d_,s,x);
227  sgnew = -d_->dot(g_->dual());
228  }
229  else {
230  sgnew = s.dot(g_->dual());
231  }
232  ls_ngrad++;
233 
234  if ( ((econd_ == CURVATURECONDITION_WOLFE)
235  && (sgnew >= c2_*sgold))
237  && (std::abs(sgnew) <= c2_*std::abs(sgold)))
239  && (c2_*sgold <= sgnew && sgnew <= -c3_*sgold))
241  && (c2_*sgold <= sgnew && sgnew <= (two*c1_ - one)*sgold)) ) {
242  curvcond = true;
243  }
244  }
245  }
246 
247  if(fnew<fmin_) {
248  fmin_ = fnew;
249  alphaMin_ = alpha;
250  }
251 
252  if (type == LINESEARCH_BACKTRACKING || type == LINESEARCH_CUBICINTERP) {
253  if (edesc_ == DESCENT_NONLINEARCG) {
254  return ((armijo && curvcond) || itcond_);
255  }
256  else {
257  return (armijo || itcond_);
258  }
259  }
260  else {
261  return ((armijo && curvcond) || itcond_);
262  }
263  }
264 
265  virtual Real getInitialAlpha(int &ls_neval, int &ls_ngrad, const Real fval, const Real gs,
266  const Vector<Real> &x, const Vector<Real> &s,
268  Real val(1), one(1), half(0.5), p1(1.e-1);
269  if (useralpha_ || usePrevAlpha_ ) {
270  val = alpha0_;
271  }
272  else {
274  Real tol = std::sqrt(ROL_EPSILON<Real>());
275  // Evaluate objective at x + s
276  updateIterate(*d_,x,s,one,con);
277  obj.update(*d_);
278  Real fnew = obj.value(*d_,tol);
279  ls_neval++;
280  // Minimize quadratic interpolate to compute new alpha
281  Real denom = (fnew - fval - gs);
282  Real alpha = ((denom > ROL_EPSILON<Real>()) ? -half*gs/denom : one);
283  val = ((alpha > p1) ? alpha : one);
284 
285  alpha0_ = val;
286  useralpha_ = true;
287  }
288  else {
289  val = one;
290  }
291  }
292  return val;
293  }
294 
295  void setNextInitialAlpha( Real alpha ) {
296  if( usePrevAlpha_ ) {
297  alpha0_ = alpha;
298  }
299  }
300 
301  void updateIterate(Vector<Real> &xnew, const Vector<Real> &x, const Vector<Real> &s, Real alpha,
302  BoundConstraint<Real> &con ) {
303 
304  xnew.set(x);
305  xnew.axpy(alpha,s);
306 
307  if ( con.isActivated() ) {
308  con.project(xnew);
309  }
310 
311  }
312 
314  return itcond_ && acceptMin_;
315  }
316 
317  bool takeNoStep() {
318  return itcond_ && !acceptMin_;
319  }
320 };
321 
322 }
323 
324 #include "ROL_LineSearchFactory.hpp"
325 
326 #endif
Provides the interface to evaluate objective functions.
void setData(Real &eps, const Vector< Real > &g)
void updateIterate(Vector< Real > &xnew, const Vector< Real > &x, const Vector< Real > &s, Real alpha, BoundConstraint< Real > &con)
virtual Real getInitialAlpha(int &ls_neval, int &ls_ngrad, const Real fval, const Real gs, const Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
Teuchos::RCP< Vector< Real > > xtst_
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Contains definitions of custom data types in ROL.
Teuchos::RCP< Vector< Real > > g_
void pruneInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=0)
Set variables to zero if they correspond to the -inactive set.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ELineSearch
Enumeration of line-search types.
Definition: ROL_Types.hpp:661
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
LineSearch(Teuchos::ParameterList &parlist)
virtual Real dot(const Vector &x) const =0
Compute where .
virtual void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=0)
Set variables to zero if they correspond to the -active set.
EDescent StringToEDescent(std::string s)
Definition: ROL_Types.hpp:402
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
bool isActivated(void)
Check if bounds are on.
Provides interface for and implements line searches.
void setMaxitUpdate(Real &alpha, Real &fnew, const Real &fold)
Teuchos::RCP< Vector< Real > > d_
Provides the interface to apply upper and lower bound constraints.
virtual ~LineSearch()
ECurvatureCondition
Enumeration of line-search curvature conditions.
Definition: ROL_Types.hpp:744
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
ECurvatureCondition StringToECurvatureCondition(std::string s)
Definition: ROL_Types.hpp:804
virtual void 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)=0
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
virtual bool status(const ELineSearch type, int &ls_neval, int &ls_ngrad, const Real alpha, const Real fold, const Real sgold, const Real fnew, const Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con)
ECurvatureCondition econd_
Teuchos::RCP< Vector< Real > > grad_
EDescent
Enumeration of descent direction types.
Definition: ROL_Types.hpp:345
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con)
void setNextInitialAlpha(Real alpha)