ROL
ROL_InteriorPointStep.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_INTERIORPOINTSTEP_H
45 #define ROL_INTERIORPOINTSTEP_H
46 
47 #include "ROL_CompositeStep.hpp"
49 #include "ROL_InteriorPoint.hpp"
51 #include "ROL_Types.hpp"
52 
53 
54 namespace ROL {
55 
56 template <class Real>
57 class InteriorPointStep : public Step<Real> {
58 
61 
63 typedef typename PV::size_type size_type;
64 
65 const static size_type OPT = 0;
66 const static size_type SLACK = 1;
67 
68 private:
69 
70  Teuchos::RCP<StatusTest<Real> > status_;
71  Teuchos::RCP<Step<Real> > step_;
72  Teuchos::RCP<IPOBJ> ipobj_;
73  Teuchos::RCP<IPCON> ipcon_;
74  Teuchos::RCP<Algorithm<Real> > algo_;
75  Teuchos::RCP<Teuchos::ParameterList> parlist_;
76 
77  // Storage
78  Teuchos::RCP<PV> x_;
79  Teuchos::RCP<Vector<Real> > g_;
80  Teuchos::RCP<Vector<Real> > l_;
81  Teuchos::RCP<Vector<Real> > c_;
82 
83  Real mu_; // Barrier parameter
84  Real mumin_; // Minimal value of barrier parameter
85  Real mumax_; // Maximal value of barrier parameter
86  Real rho_; // Barrier parameter reduction factor
87  int maxit_; // Maximum number of interior point subproblem solves
88 
89  // For the subproblem
90  Real gtol_; // Status test gradient tolerance
91  Real ctol_; // Status test constraint tolerance
92  Real stol_; // Status test step tolerance
93  int subproblemIter_; // Status test maximum number of iterations
94 
95  int verbosity_; // Adjust level of detail in printing step information
96 
97 public:
98 
100  using Step<Real>::compute;
101  using Step<Real>::update;
102 
104 
105  InteriorPointStep(Teuchos::ParameterList &parlist) :
106  Step<Real>(),
107  status_(Teuchos::null),
108  step_(Teuchos::null),
109  ipobj_(Teuchos::null),
110  ipcon_(Teuchos::null),
111  algo_(Teuchos::null),
112  x_(Teuchos::null),
113  g_(Teuchos::null),
114  l_(Teuchos::null),
115  c_(Teuchos::null) {
116 
117  using Teuchos::ParameterList;
118 
119  verbosity_ = parlist.sublist("General").get("Print Verbosity",0);
120 
121  // List of general Interior Point parameters
122  ParameterList& iplist = parlist.sublist("Step").sublist("Interior Point");
123 
124  mu_ = iplist.get("Initial Barrier Penalty",1.0);
125  mumin_ = iplist.get("Minimum Barrier Penalty",1.e-4);
126  mumax_ = iplist.get("Maximum Barrier Penalty",1e8);
127  rho_ = iplist.get("Barrier Penalty Reduction Factor",0.5);
128  subproblemIter_ = iplist.get("Subproblem Iteration Limit",10);
129 
130 
131  // List of Status Test parameters
132  ParameterList& stlist = parlist.sublist("Status Test");
133 
134  gtol_ = stlist.get("Gradient Tolerance", 1.e-8);
135  ctol_ = stlist.get("Constraint Tolerance", 1.e-8);
136  stol_ = stlist.get("Step Tolerance", 1.e-8);
137  maxit_ = stlist.get("Iteration Limit", 100);
138 
139  parlist_ = Teuchos::rcp(&parlist, false);
140 
141  }
142 
145  void initialize( Vector<Real> &x, const Vector<Real> &g,
146  Vector<Real> &l, const Vector<Real> &c,
148  AlgorithmState<Real> &algo_state ) {
149 
150  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
151  state->descentVec = x.clone();
152  state->gradientVec = g.clone();
153  state->constraintVec = c.clone();
154 
155  // Initialize storage
156  x_ = Teuchos::rcp_static_cast<PV>(x.clone());
157  g_ = g.clone();
158  l_ = l.clone();
159  c_ = c.clone();
160 
161  x_->set(x);
162 
163  ipobj_ = Teuchos::rcp(&Teuchos::dyn_cast<IPOBJ>(obj),false);
164  ipcon_ = Teuchos::rcp(&Teuchos::dyn_cast<IPCON>(con),false);
165 
166  // Set initial penalty
167  ipobj_->updatePenalty(mu_);
168 
169  algo_state.nfval = 0;
170  algo_state.ncval = 0;
171  algo_state.ngrad = 0;
172 
173  Real zerotol = 0.0;
174  obj.update(*x_,true,algo_state.iter);
175  algo_state.value = obj.value(*x_,zerotol);
176 
177  obj.gradient(*g_,*x_,zerotol);
178  algo_state.gnorm = g_->norm();
179 
180  con.value(*c_,*x_,zerotol);
181  algo_state.cnorm = c_->norm();
182 
183  algo_state.nfval += ipobj_->getNumberFunctionEvaluations();
184  algo_state.ngrad += ipobj_->getNumberGradientEvaluations();
185  algo_state.ncval += ipcon_->getNumberConstraintEvaluations();
186 
187  }
188 
189 
190 
193  AlgorithmState<Real> &algo_state ) {
194  initialize(x,g,l,c,obj,con,algo_state);
195  }
196 
197 
198 
199 
202  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
204  AlgorithmState<Real> &algo_state ) {
205 
206  // Create the algorithm
207  algo_ = Teuchos::rcp( new Algorithm<Real>("Composite Step",*parlist_,false) );
208 
209  x_->set(x);
210 
211  // Run the algorithm
212  algo_->run(*x_,*g_,*l_,*c_,*ipobj_,*ipcon_,false);
213 
214  s.set(*x_); s.axpy(-1.0,x);
215 
216  // Get number of iterations from the subproblem solve
217  subproblemIter_ = (algo_->getState())->iter;
218 
219  }
220 
221  virtual void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
224  AlgorithmState<Real> &algo_state ) {
225  compute(s,x,l,obj,con,algo_state);
226  }
227 
228 
229 
233  EqualityConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
234 
235  // If we can change the barrier parameter, do so
236  if( (rho_< 1.0 && mu_ > mumin_) || (rho_ > 1.0 && mu_ < mumax_) ) {
237  mu_ *= rho_;
238  ipobj_->updatePenalty(mu_);
239  }
240 
241  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
242 
243  // Update optimization vector
244  x.plus(s);
245 
246  algo_state.iterateVec->set(x);
247  state->descentVec->set(s);
248  algo_state.snorm = s.norm();
249  algo_state.iter++;
250 
251  Real zerotol = 0.0;
252 
253  algo_state.value = ipobj_->value(x,zerotol);
254  algo_state.value = ipobj_->getObjectiveValue();
255 
256  ipcon_->value(*c_,x,zerotol);
257  state->constraintVec->set(*c_);
258 
259  ipobj_->gradient(*g_,x,zerotol);
260  state->gradientVec->set(*g_);
261 
262  ipcon_->applyAdjointJacobian(*g_,*l_,x,zerotol);
263  state->gradientVec->plus(*g_);
264 
265  x_->set(x);
266  x_->axpy(-1.0,state->gradientVec->dual());
267 
268  Elementwise::ThresholdUpper<Real> threshold(0.0);
269 
270  //PartitionedVector<Real> &xpv = Teuchos::dyn_cast<PartitionedVector<Real> >(*x_);
271 
272  Teuchos::RCP<Vector<Real> > slack = x_->get(SLACK);
273 
274  slack->applyUnary(threshold);
275 
276  x_->axpy(-1.0,x);
277 
278  algo_state.gnorm = x_->norm();
279  algo_state.cnorm = state->constraintVec->norm();
280  algo_state.snorm = s.norm();
281 
282  algo_state.nfval += ipobj_->getNumberFunctionEvaluations();
283  algo_state.ngrad += ipobj_->getNumberGradientEvaluations();
284  algo_state.ncval += ipcon_->getNumberConstraintEvaluations();
285 
286  }
287 
291  AlgorithmState<Real> &algo_state ) {
292  update(x,l,s,obj,con,algo_state);
293  }
294 
295 
296 
300  virtual void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj,
301  BoundConstraint<Real> &con,
302  AlgorithmState<Real> &algo_state ) {}
303 
307  virtual void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj,
309  AlgorithmState<Real> &algo_state ) {}
310 
313  std::string printHeader( void ) const {
314  std::stringstream hist;
315 
316  if( verbosity_ > 0 ) {
317 
318  hist << std::string(116,'-') << "\n";
319  hist << "Interior Point status output definitions\n\n";
320 
321  hist << " IPiter - Number of interior point steps taken\n";
322  hist << " CSiter - Number of Composite Steps taken in each subproblem\n";
323  hist << " penalty - Penalty parameter multiplying the barrier objective\n";
324  hist << " fval - Number of objective evaluations\n";
325  hist << " cnorm - Norm of the composite constraint\n";
326  hist << " gLnorm - Norm of the Lagrangian's gradient\n";
327  hist << " snorm - Norm of step (update to optimzation and slack vector)\n";
328  hist << " #fval - Number of objective function evaluations\n";
329  hist << " #grad - Number of gradient evaluations\n";
330  hist << " #cval - Number of composite constraint evaluations\n";
331  hist << std::string(116,'-') << "\n";
332 
333 
334  }
335 
336  hist << " ";
337  hist << std::setw(9) << std::left << "IPiter";
338  hist << std::setw(9) << std::left << "CSiter";
339  hist << std::setw(15) << std::left << "penalty";
340  hist << std::setw(15) << std::left << "fval";
341  hist << std::setw(15) << std::left << "cnorm";
342  hist << std::setw(15) << std::left << "gLnorm";
343  hist << std::setw(15) << std::left << "snorm";
344  hist << std::setw(8) << std::left << "#fval";
345  hist << std::setw(8) << std::left << "#grad";
346  hist << std::setw(8) << std::left << "#cval";
347 
348  hist << "\n";
349  return hist.str();
350  }
351 
354  std::string printName( void ) const {
355  std::stringstream hist;
356  hist << "\n" << "Composite Step Interior Point Solver\n";
357  return hist.str();
358  }
359 
362  std::string print( AlgorithmState<Real> &algo_state, bool pHeader = false ) const {
363  std::stringstream hist;
364  hist << std::scientific << std::setprecision(6);
365  if ( algo_state.iter == 0 ) {
366  hist << printName();
367  }
368  if ( pHeader ) {
369  hist << printHeader();
370  }
371  if ( algo_state.iter == 0 ) {
372  hist << " ";
373  hist << std::setw(9) << std::left << algo_state.iter;
374  hist << std::setw(9) << std::left << subproblemIter_;
375  hist << std::setw(15) << std::left << mu_;
376  hist << std::setw(15) << std::left << algo_state.value;
377  hist << std::setw(15) << std::left << algo_state.cnorm;
378  hist << std::setw(15) << std::left << algo_state.gnorm;
379  hist << "\n";
380  }
381  else {
382  hist << " ";
383  hist << std::setw(9) << std::left << algo_state.iter;
384  hist << std::setw(9) << std::left << subproblemIter_;
385  hist << std::setw(15) << std::left << mu_;
386  hist << std::setw(15) << std::left << algo_state.value;
387  hist << std::setw(15) << std::left << algo_state.cnorm;
388  hist << std::setw(15) << std::left << algo_state.gnorm;
389  hist << std::setw(15) << std::left << algo_state.snorm;
390 // hist << std::scientific << std::setprecision(6);
391  hist << std::setw(8) << std::left << algo_state.nfval;
392  hist << std::setw(8) << std::left << algo_state.ngrad;
393  hist << std::setw(8) << std::left << algo_state.ncval;
394  hist << "\n";
395  }
396 
397  return hist.str();
398  }
399 
400 
401 
402 
403 
404 }; // class InteriorPointStep
405 
406 } // namespace ROL
407 
408 #endif // ROL_INTERIORPOINTSTEP_H
Provides the interface to evaluate objective functions.
std::string printName(void) const
Print step name.
Teuchos::RCP< Teuchos::ParameterList > parlist_
virtual void plus(const Vector &x)=0
Compute , where .
InteriorPoint::PenalizedObjective< Real > IPOBJ
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, EqualityConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with equality constraint.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
InteriorPoint::CompositeConstraint< Real > IPCON
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
static const size_type SLACK
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:69
Defines the linear algebra of vector space on a generic partitioned vector.
Has both inequality and equality constraints. Treat inequality constraint as equality with slack vari...
Teuchos::RCP< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:74
Contains definitions of custom data types in ROL.
Teuchos::RCP< Step< Real > > step_
Teuchos::RCP< IPOBJ > ipobj_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
virtual void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, for bound constraints; here only to satisfy the interface requirements, does nothing, needs refactoring.
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, EqualityConstraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Initialize step with equality constraint.
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:91
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, EqualityConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step (equality constraints).
Defines the equality constraint operator interface.
Teuchos::RCP< Vector< Real > > g_
PartitionedVector< Real > PV
Provides an interface to run optimization algorithms.
virtual void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step for bound constraints; here only to satisfy the interface requirements, does nothing, needs refactoring.
Provides the interface to apply upper and lower bound constraints.
Teuchos::RCP< Algorithm< Real > > algo_
Teuchos::RCP< Vector< Real > > c_
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, EqualityConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful (equality constraints).
std::vector< PV >::size_type size_type
static const size_type OPT
virtual void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, EqualityConstraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Compute step (equality constraints).
Teuchos::RCP< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:105
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
virtual Real norm() const =0
Returns where .
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, EqualityConstraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Update step, if successful (equality constraints).
Teuchos::RCP< Vector< Real > > l_
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Teuchos::RCP< IPCON > ipcon_
std::string printHeader(void) const
Print iterate header.
Teuchos::RCP< StatusTest< Real > > status_
InteriorPointStep(Teuchos::ParameterList &parlist)