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