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 eps_; // Minimal value of barrier parameter
85  Real rho_; // Barrier parameter reduction factor
86  int maxit_; // Maximum number of interior point subproblem solves
87 
88  // For the subproblem
89  Real gtol_; // Status test gradient tolerance
90  Real ctol_; // Status test constraint tolerance
91  Real stol_; // Status test step tolerance
92  int subproblemIter_; // Status test maximum number of iterations
93 
94 public:
95 
97 
98  InteriorPointStep(Teuchos::ParameterList &parlist) :
99  Step<Real>(),
100  status_(Teuchos::null),
101  step_(Teuchos::null),
102  ipobj_(Teuchos::null),
103  ipcon_(Teuchos::null),
104  algo_(Teuchos::null),
105  x_(Teuchos::null),
106  g_(Teuchos::null),
107  l_(Teuchos::null),
108  c_(Teuchos::null) {
109 
110  using Teuchos::ParameterList;
111 
112  // List of general Interior Point parameters
113  ParameterList& iplist = parlist.sublist("Step").sublist("Interior Point");
114 
115  mu_ = iplist.get("Initial Barrier Penalty",1.0);
116  eps_ = iplist.get("Minimum Barrier Penalty",1.e-4);
117  rho_ = iplist.get("Barrier Penalty Reduction Factor",0.5);
118  subproblemIter_ = iplist.get("Subproblem Iteration Limit",10);
119 
120 
121  // List of Status Test parameters
122  ParameterList& stlist = parlist.sublist("Status Test");
123 
124  gtol_ = stlist.get("Gradient Tolerance", 1.e-8);
125  ctol_ = stlist.get("Constraint Tolerance", 1.e-8);
126  stol_ = stlist.get("Step Tolerance", 1.e-8);
127  maxit_ = stlist.get("Iteration Limit", 100);
128 
129  parlist_ = Teuchos::rcp(&parlist, false);
130 
131  }
132 
135  void initialize( Vector<Real> &x, const Vector<Real> &g,
136  Vector<Real> &l, const Vector<Real> &c,
138  AlgorithmState<Real> &algo_state ) {
139 
140  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
141  state->descentVec = x.clone();
142  state->gradientVec = g.clone();
143  state->constraintVec = c.clone();
144 
145  // Initialize storage
146  x_ = Teuchos::rcp_static_cast<PV>(x.clone());
147  g_ = g.clone();
148  l_ = l.clone();
149  c_ = c.clone();
150 
151  x_->set(x);
152 
153  ipobj_ = Teuchos::rcp(&Teuchos::dyn_cast<IPOBJ>(obj),false);
154  ipcon_ = Teuchos::rcp(&Teuchos::dyn_cast<IPCON>(con),false);
155 
156  // Set initial penalty
157  ipobj_->updatePenalty(mu_);
158 
159  algo_state.nfval = 0;
160  algo_state.ncval = 0;
161  algo_state.ngrad = 0;
162 
163  Real zerotol = 0.0;
164  obj.update(*x_,true,algo_state.iter);
165  algo_state.value = obj.value(*x_,zerotol);
166 
167  obj.gradient(*g_,*x_,zerotol);
168  algo_state.gnorm = g_->norm();
169 
170  con.value(*c_,*x_,zerotol);
171  algo_state.cnorm = c_->norm();
172 
173  algo_state.nfval += ipobj_->getNumberFunctionEvaluations();
174  algo_state.ngrad += ipobj_->getNumberGradientEvaluations();
175  algo_state.ncval += ipcon_->getNumberConstraintEvaluations();
176 
177  }
178 
179 
180 
183  AlgorithmState<Real> &algo_state ) {
184  initialize(x,g,l,c,obj,con,algo_state);
185  }
186 
187 
188 
189 
192  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
194  AlgorithmState<Real> &algo_state ) {
195 
196  // Create the algorithm
197  algo_ = Teuchos::rcp( new Algorithm<Real>("Composite Step",*parlist_,false) );
198 
199  x_->set(x);
200 
201  // Run the algorithm
202  algo_->run(*x_,*g_,*l_,*c_,*ipobj_,*ipcon_,false);
203 
204  s.set(*x_); s.axpy(-1.0,x);
205 
206  // Get number of iterations from the subproblem solve
207  subproblemIter_ = (algo_->getState())->iter;
208 
209  }
210 
211  virtual void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
214  AlgorithmState<Real> &algo_state ) {
215  compute(s,x,l,obj,con,algo_state);
216  }
217 
218 
219 
223  EqualityConstraint<Real> &con, AlgorithmState<Real> &algo_state ) {
224 
225  // If we can reduce the barrier parameter, do so
226  if(mu_ > eps_) {
227  mu_ *= rho_;
228  ipobj_->updatePenalty(mu_);
229  }
230 
231  Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
232 
233  // Update optimization vector
234  x.plus(s);
235 
236  algo_state.iterateVec->set(x);
237  state->descentVec->set(s);
238  algo_state.snorm = s.norm();
239  algo_state.iter++;
240 
241  Real zerotol = 0.0;
242 
243  algo_state.value = ipobj_->value(x,zerotol);
244  algo_state.value = ipobj_->getObjectiveValue();
245 
246  ipcon_->value(*c_,x,zerotol);
247  state->constraintVec->set(*c_);
248 
249  ipobj_->gradient(*g_,x,zerotol);
250  state->gradientVec->set(*g_);
251 
252  ipcon_->applyAdjointJacobian(*g_,*l_,x,zerotol);
253  state->gradientVec->plus(*g_);
254 
255  x_->set(x);
256  x_->axpy(-1.0,state->gradientVec->dual());
257 
258  Elementwise::ThresholdUpper<Real> threshold(0.0);
259 
260  //PartitionedVector<Real> &xpv = Teuchos::dyn_cast<PartitionedVector<Real> >(*x_);
261 
262  Teuchos::RCP<Vector<Real> > slack = x_->get(SLACK);
263 
264  slack->applyUnary(threshold);
265 
266  x_->axpy(-1.0,x);
267 
268  algo_state.gnorm = x_->norm();
269  algo_state.cnorm = state->constraintVec->norm();
270  algo_state.snorm = s.norm();
271 
272  algo_state.nfval += ipobj_->getNumberFunctionEvaluations();
273  algo_state.ngrad += ipobj_->getNumberGradientEvaluations();
274  algo_state.ncval += ipcon_->getNumberConstraintEvaluations();
275 
276  }
277 
281  AlgorithmState<Real> &algo_state ) {
282  update(x,l,s,obj,con,algo_state);
283  }
284 
285 
286 
290  virtual void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj,
291  BoundConstraint<Real> &con,
292  AlgorithmState<Real> &algo_state ) {}
293 
297  virtual void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj,
299  AlgorithmState<Real> &algo_state ) {}
300 
303  std::string printHeader( void ) const {
304  std::stringstream hist;
305  hist << " ";
306  hist << std::setw(9) << std::left << "IPiter";
307  hist << std::setw(9) << std::left << "CSiter";
308  hist << std::setw(15) << std::left << "penalty";
309  hist << std::setw(15) << std::left << "fval";
310  hist << std::setw(15) << std::left << "cnorm";
311  hist << std::setw(15) << std::left << "gLnorm";
312  hist << std::setw(15) << std::left << "snorm";
313  hist << std::setw(8) << std::left << "#fval";
314  hist << std::setw(8) << std::left << "#grad";
315  hist << std::setw(8) << std::left << "#cval";
316 
317  hist << "\n";
318  return hist.str();
319  }
320 
323  std::string printName( void ) const {
324  std::stringstream hist;
325  hist << "\n" << "Composite Step Interior Point Solver\n";
326  return hist.str();
327  }
328 
331  std::string print( AlgorithmState<Real> &algo_state, bool pHeader = false ) const {
332  std::stringstream hist;
333  hist << std::scientific << std::setprecision(6);
334  if ( algo_state.iter == 0 ) {
335  hist << printName();
336  }
337  if ( pHeader ) {
338  hist << printHeader();
339  }
340  if ( algo_state.iter == 0 ) {
341  hist << " ";
342  hist << std::setw(9) << std::left << algo_state.iter;
343  hist << std::setw(9) << std::left << subproblemIter_;
344  hist << std::setw(15) << std::left << mu_;
345  hist << std::setw(15) << std::left << algo_state.value;
346  hist << std::setw(15) << std::left << algo_state.cnorm;
347  hist << std::setw(15) << std::left << algo_state.gnorm;
348  hist << "\n";
349  }
350  else {
351  hist << " ";
352  hist << std::setw(9) << std::left << algo_state.iter;
353  hist << std::setw(9) << std::left << subproblemIter_;
354  hist << std::setw(15) << std::left << mu_;
355  hist << std::setw(15) << std::left << algo_state.value;
356  hist << std::setw(15) << std::left << algo_state.cnorm;
357  hist << std::setw(15) << std::left << algo_state.gnorm;
358  hist << std::setw(15) << std::left << algo_state.snorm;
359 // hist << std::scientific << std::setprecision(6);
360  hist << std::setw(8) << std::left << algo_state.nfval;
361  hist << std::setw(8) << std::left << algo_state.ngrad;
362  hist << std::setw(8) << std::left << algo_state.ncval;
363  hist << "\n";
364  }
365 
366  return hist.str();
367  }
368 
369 
370 
371 
372 
373 }; // class InteriorPointStep
374 
375 } // namespace ROL
376 
377 #endif // ROL_INTERIORPOINTSTEP_H
Provides the interface to evaluate objective functions.
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
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:67
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:72
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:77
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.
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
Provides the interface to apply upper and lower bound constraints.
std::string printName(void) const
Print step name.
Teuchos::RCP< Algorithm< Real > > algo_
std::string printHeader(void) const
Print iterate header.
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:91
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_
Teuchos::RCP< StatusTest< Real > > status_
InteriorPointStep(Teuchos::ParameterList &parlist)