ROL
ROL_InteriorPointObjective.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_INTERIORPOINTOBJECTIVE_H
11 #define ROL_INTERIORPOINTOBJECTIVE_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_BoundConstraint.hpp"
15 #include "ROL_ParameterList.hpp"
16 #include "ROL_ScalarController.hpp"
17 
27 namespace ROL {
28 
29 template<class Real>
30 class InteriorPointObjective : public Objective<Real> {
31 
32  typedef Elementwise::ValueSet<Real> ValueSet;
33 
34 private:
35 
36  const Ptr<Objective<Real>> obj_;
37  const Ptr<BoundConstraint<Real>> bnd_;
38  const Ptr<const Vector<Real>> lo_;
39  const Ptr<const Vector<Real>> up_;
40 
41  Ptr<Vector<Real>> maskL_; // Elements are 1 when xl>-INF, zero for xl =-INF
42  Ptr<Vector<Real>> maskU_; // Elements are 1 when xu< INF, zero for xu = INF
43  Ptr<Vector<Real>> maskL0_; // Elements are 1 when xl>-INF and xu = INF, zero for xl =-INF
44  Ptr<Vector<Real>> maskU0_; // Elements are 1 when xu< INF and XL =-INF, zero for xu = INF
45  Ptr<Vector<Real>> pwa_; // Scratch vector
46 
47  bool useLinearDamping_; // Add linear damping terms to the penalized objective
48  // to prevent the problems such as when the log barrier
49  // contribution is unbounded below on the feasible set
50  Real kappaD_; // Linear damping coefficient
51  Real mu_; // Penalty parameter
52 
53  Ptr<ScalarController<Real,int>> fval_;
54  Ptr<VectorController<Real,int>> gradient_;
55 
56  int nfval_;
57  int ngrad_;
58 
59  // x <- f(x) = { log(x) if x > 0
60  // { -inf if x <= 0
61  class ModifiedLogarithm : public Elementwise::UnaryFunction<Real> {
62  public:
63  Real apply( const Real &x ) const {
64  const Real zero(0), NINF(ROL_NINF<Real>());
65  return (x>zero) ? std::log(x) : NINF;
66  //return std::log(x);
67  }
68  }; // class ModifiedLogarithm
69 
70  // x <- f(x) = { 1/x if x > 0
71  // { 0 if x <= 0
72  class ModifiedReciprocal : public Elementwise::UnaryFunction<Real> {
73  public:
74  Real apply( const Real &x ) const {
75  const Real zero(0), one(1);
76  return (x>zero) ? one/x : zero;
77  //return one/x;
78  }
79 
80  }; // class ModifiedReciprocal
81 
82  // x <- f(x,y) = { y/x if x > 0
83  // { 0 if x <= 0
84  class ModifiedDivide : public Elementwise::BinaryFunction<Real> {
85  public:
86  Real apply( const Real &x, const Real &y ) const {
87  const Real zero(0);
88  return (x>zero) ? y/x : zero;
89  //return y/x;
90  }
91  }; // class ModifiedDivide
92 
93  // x <- f(x,y) = { x if y != 0, complement == false
94  // { 0 if y == 0, complement == false
95  // { 0 if y != 0, complement == true
96  // { x if y == 0, complement == true
97  class Mask : public Elementwise::BinaryFunction<Real> {
98  private:
100  public:
101  Mask( bool complement ) : complement_(complement) {}
102  Real apply( const Real &x, const Real &y ) const {
103  const Real zero(0);
104  return ( complement_ ^ (y != zero) ) ? zero : x;
105  }
106  }; // class Mask
107 
108  void initialize(const Vector<Real> &x, const Vector<Real> &g) {
109  const Real zero(0), one(1);
110 
111  fval_ = makePtr<ScalarController<Real,int>>();
112  gradient_ = makePtr<VectorController<Real,int>>();
113 
114  // Determine the index sets where the
115  ValueSet isBoundedBelow( ROL_NINF<Real>(), ValueSet::GREATER_THAN, one, zero );
116  ValueSet isBoundedAbove( ROL_INF<Real>(), ValueSet::LESS_THAN, one, zero );
117 
118  maskL_ = x.clone(); maskL_->applyBinary(isBoundedBelow,*lo_);
119  maskU_ = x.clone(); maskU_->applyBinary(isBoundedAbove,*up_);
120 
121  pwa_ = x.clone();
122 
123  if( useLinearDamping_ ) {
124  maskL0_ = x.clone();
125  maskL0_->set(*maskL_); // c_i = { 1 if l_i > NINF
126  // { 0 otherwise
127  maskL0_->applyBinary(Mask(true),*maskU_); // c_i = { 1 if l_i > NINF and u_i = INF
128  // { 0 otherwise
129  maskU0_ = x.clone();
130  maskU0_->set(*maskU_); // c_i = { 1 if u_i < INF
131  // { 0 otherwise
132  maskU0_->applyBinary(Mask(true),*maskL_); // c_i = { 1 if u_i < INF and l_i = NINF
133  // { 0 otherwise
134  }
135  }
136 
137 public:
138 
140  const Ptr<BoundConstraint<Real>> &bnd,
141  const Vector<Real> &x,
142  const Vector<Real> &g,
143  const bool useLinearDamping,
144  const Real kappaD,
145  const Real mu )
146  : obj_(obj), bnd_(bnd), lo_(bnd->getLowerBound()), up_(bnd->getUpperBound()),
147  useLinearDamping_(useLinearDamping), kappaD_(kappaD), mu_(mu),
148  nfval_(0), ngrad_(0) {
149  initialize(x,g);
150  }
151 
153  const Ptr<BoundConstraint<Real>> &bnd,
154  const Vector<Real> &x,
155  const Vector<Real> &g,
156  ParameterList &parlist )
157  : obj_(obj), bnd_(bnd), lo_(bnd->getLowerBound()), up_(bnd->getUpperBound()),
158  nfval_(0), ngrad_(0) {
159  ParameterList &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
160  ParameterList &lblist = iplist.sublist("Barrier Objective");
161 
162  useLinearDamping_ = lblist.get("Use Linear Damping", true);
163  kappaD_ = lblist.get("Linear Damping Coefficient", 1.e-4);
164  mu_ = lblist.get("Initial Barrier Parameter", 0.1);
165 
166  initialize(x,g);
167  }
168 
169  Real getObjectiveValue(const Vector<Real> &x, Real &tol) {
170  int key(0);
171  Real val(0);
172  bool isComputed = fval_->get(val,key);
173  if (!isComputed) {
174  val = obj_->value(x,tol); nfval_++;
175  fval_->set(val,key);
176  }
177  return val;
178  }
179 
180  const Ptr<const Vector<Real>> getObjectiveGradient(const Vector<Real> &x, Real &tol) {
181  int key(0);
182  if (!gradient_->isComputed(key)) {
183  if (gradient_->isNull(key)) gradient_->allocate(x.dual(),key);
184  obj_->gradient(*gradient_->set(key),x,tol); ngrad_++;
185  }
186  return gradient_->get(key);
187  }
188 
190  return nfval_;
191  }
192 
194  return ngrad_;
195  }
196 
197  void updatePenalty(const Real mu) {
198  mu_ = mu;
199  }
200 
201  void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
202  obj_->update(x,type,iter);
203  fval_->objectiveUpdate(type);
204  gradient_->objectiveUpdate(type);
205  }
206 
207  Real value( const Vector<Real> &x, Real &tol ) {
208  const Real zero(0), one(1);
209  Real linearTerm = zero;
210  // Compute the unpenalized objective value
211  Real fval = getObjectiveValue(x,tol);
212  // Compute log barrier
213  ModifiedLogarithm mlog;
214  Elementwise::ReductionSum<Real> sum;
215  Elementwise::Multiply<Real> mult;
216 
217  pwa_->set(x); // pwa = x
218  pwa_->axpy(-one,*lo_); // pwa = x-l
219  if( useLinearDamping_ ) {
220  // Penalizes large positive x_i when only a lower bound exists
221  linearTerm += maskL0_->dot(*pwa_);
222  }
223  pwa_->applyUnary(mlog); // pwa = mlog(x-l)
224  Real aval = pwa_->dot(*maskL_);
225 
226  pwa_->set(*up_); // pwa = u
227  pwa_->axpy(-one,x); // pwa = u-x
228  if( useLinearDamping_ ) {
229  // Penalizes large negative x_i when only an upper bound exists
230  linearTerm += maskU0_->dot(*pwa_);
231  }
232  pwa_->applyUnary(mlog); // pwa = mlog(u-x)
233  Real bval = pwa_->dot(*maskU_);
234 
235  fval -= mu_*(aval+bval);
236  fval += kappaD_*mu_*linearTerm;
237  return fval;
238  }
239 
240  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
241  const Real one(1);
242  // Compute gradient of objective function
243  g.set(*getObjectiveGradient(x,tol));
244 
245  // Add gradient of the log barrier penalty
246  ModifiedReciprocal mrec;
247 
248  pwa_->set(x); // pwa = x
249  pwa_->axpy(-one,*lo_); // pwa = x-l
250  pwa_->applyUnary(mrec); // pwa_i = 1/(x_i-l_i) for i s.t. x_i > l_i, 0 otherwise
251  pwa_->applyBinary(Mask(true),*maskL_); // zero elements where l = NINF
252  g.axpy(-mu_,pwa_->dual());
253  if( useLinearDamping_ ) {
254  g.axpy(-mu_*kappaD_,maskL0_->dual());
255  }
256 
257  pwa_->set(*up_); // pwa = u
258  pwa_->axpy(-one,x); // pwa = u-x
259  pwa_->applyUnary(mrec); // pwa_i = 1/(u_i-x_i) for i s.t. u_i > x_i, 0 otherwise
260  pwa_->applyBinary(Mask(true),*maskU_); // zero elements where u = INF
261  g.axpy( mu_,pwa_->dual());
262  if( useLinearDamping_ ) {
263  g.axpy( mu_*kappaD_,maskU0_->dual());
264  }
265  }
266 
267  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
268  const Real one(1), two(2);
269  // Evaluate objective hessian
270  obj_->hessVec(hv,v,x,tol);
271 
272  // Evaluate log barrier hessian
273  ModifiedReciprocal mrec;
274  Elementwise::Multiply<Real> mult;
275  Elementwise::Power<Real> square(two);
276 
277  pwa_->set(x); // pwa = x
278  pwa_->axpy(-one,*lo_); // pwa = x-l
279  pwa_->applyUnary(mrec); // pwa_i = 1/(x_i-l_i) for i s.t. x_i > l_i, 0 otherwise
280  pwa_->applyBinary(Mask(true),*maskL_); // zero elements where l = NINF
281  pwa_->applyUnary(square); // pwa_i = { (x_i-l_i)^(-2) if l_i > NINF
282  // { 0 if l_i = NINF
283  pwa_->applyBinary(mult,v);
284  hv.axpy(mu_,pwa_->dual());
285 
286  pwa_->set(*up_); // pwa = u
287  pwa_->axpy(-one,x); // pwa = u-x
288  pwa_->applyUnary(mrec); // pwa_i = 1/(u_i-x_i) for i s.t. u_i > x_i, 0 otherwise
289  pwa_->applyBinary(Mask(true),*maskU_); // zero elements where u = INF
290  pwa_->applyUnary(square); // pwa_i = { (u_i-x_i)^(-2) if u_i < INF
291  // { 0 if u_i = INF
292  pwa_->applyBinary(mult,v);
293  hv.axpy(mu_,pwa_->dual());
294  }
295 
296 }; // class InteriorPointObjective
297 
298 }
299 
300 #endif // ROL_INTERIORPOINTOBJECTIVE_H
Provides the interface to evaluate objective functions.
InteriorPointObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, ParameterList &parlist)
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
const Ptr< const Vector< Real > > lo_
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Ptr< ScalarController< Real, int > > fval_
Real apply(const Real &x, const Real &y) const
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
const Ptr< Objective< Real > > obj_
Provides the interface to apply upper and lower bound constraints.
Elementwise::ValueSet< Real > ValueSet
const Ptr< BoundConstraint< Real > > bnd_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
const Ptr< const Vector< Real > > getObjectiveGradient(const Vector< Real > &x, Real &tol)
const Ptr< const Vector< Real > > up_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
Ptr< VectorController< Real, int > > gradient_
InteriorPointObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, const bool useLinearDamping, const Real kappaD, const Real mu)
Real apply(const Real &x, const Real &y) const
Real getObjectiveValue(const Vector< Real > &x, Real &tol)