ROL
ROL_InteriorPointPenalty.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_INTERIORPOINTPENALTY_H
45#define ROL_INTERIORPOINTPENALTY_H
46
47#include "ROL_Objective.hpp"
49#include "ROL_ParameterList.hpp"
50
60namespace ROL {
61
62template<class Real>
63class InteriorPointPenalty : public Objective<Real> {
64
65 typedef Vector<Real> V;
68
69 typedef Elementwise::ValueSet<Real> ValueSet;
70
71private:
72
73 const ROL::Ptr<OBJ> obj_;
74 const ROL::Ptr<BND> bnd_;
75 const ROL::Ptr<const V> lo_;
76 const ROL::Ptr<const V> up_;
77
78 ROL::Ptr<V> g_; // Gradient of penalized objective
79
80 ROL::Ptr<V> maskL_; // Elements are 1 when xl>-INF, zero for xl = -INF
81 ROL::Ptr<V> maskU_; // Elements are 1 when xu< INF, zero for xu = INF
82
83 ROL::Ptr<V> a_; // Scratch vector
84 ROL::Ptr<V> b_; // Scratch vector
85 ROL::Ptr<V> c_; // Scratch vector
86
87 bool useLinearDamping_; // Add linear damping terms to the penalized objective
88 // to prevent the problems such as when the log barrier
89 // contribution is unbounded below on the feasible set
90
91 Real mu_; // Penalty parameter
92 Real kappaD_; // Linear damping coefficient
93 Real fval_; // Unpenalized objective value
94
95 int nfval_;
96 int ngval_;
97
98
99
100 // x <- f(x) = { log(x) if x > 0
101 // { 0 if x <= 0
102 class ModifiedLogarithm : public Elementwise::UnaryFunction<Real> {
103 public:
104 Real apply( const Real &x ) const {
105 return (x>0) ? std::log(x) : Real(0.0);
106 }
107 }; // class ModifiedLogarithm
108
109 // x <- f(x) = { 1/x if x > 0
110 // { 0 if x <= 0
111 class ModifiedReciprocal : public Elementwise::UnaryFunction<Real> {
112 public:
113 Real apply( const Real &x ) const {
114 return (x>0) ? 1.0/x : Real(0.0);
115 }
116
117 }; // class ModifiedReciprocal
118
119 // x <- f(x,y) = { y/x if x > 0
120 // { 0 if x <= 0
121 class ModifiedDivide : public Elementwise::BinaryFunction<Real> {
122 public:
123 Real apply( const Real &x, const Real &y ) const {
124 return (x>0) ? y/x : Real(0.0);
125 }
126 }; // class ModifiedDivide
127
128 // x <- f(x,y) = { x if y != 0, complement == false
129 // { 0 if y == 0, complement == false
130 // { 0 if y != 0, complement == true
131 // { x if y == 0, complement == true
132 class Mask : public Elementwise::BinaryFunction<Real> {
133 private:
135 public:
136 Mask( bool complement ) : complement_(complement) {}
137 Real apply( const Real &x, const Real &y ) const {
138 return ( complement_ ^ (y !=0) ) ? 0 : x;
139 }
140 }; // class Mask
141
142
143public:
144
146
147 InteriorPointPenalty( const ROL::Ptr<Objective<Real> > &obj,
148 const ROL::Ptr<BoundConstraint<Real> > &con,
149 ROL::ParameterList &parlist ) :
150 obj_(obj), bnd_(con), lo_( con->getLowerBound() ), up_( con->getUpperBound() ) {
151
152 Real one(1);
153 Real zero(0);
154
155 // Determine the index sets where the
156 ValueSet isBoundedBelow( ROL_NINF<Real>(), ValueSet::GREATER_THAN, one, zero );
157 ValueSet isBoundedAbove( ROL_INF<Real>(), ValueSet::LESS_THAN, one, zero );
158
159 maskL_ = lo_->clone();
160 maskU_ = up_->clone();
161
162 maskL_->applyBinary(isBoundedBelow,*lo_);
163 maskU_->applyBinary(isBoundedAbove,*up_);
164
165 ROL::ParameterList &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
166 ROL::ParameterList &lblist = iplist.sublist("Barrier Objective");
167
168 useLinearDamping_ = lblist.get("Use Linear Damping",true);
169 kappaD_ = lblist.get("Linear Damping Coefficient",1.e-4);
170 mu_ = lblist.get("Initial Barrier Parameter",0.1);
171
172
173 a_ = lo_->clone();
174 b_ = up_->clone();
175 g_ = lo_->dual().clone();
176
177 if( useLinearDamping_ ) {
178 c_ = lo_->clone();
179 }
180 }
181
182 Real getObjectiveValue(void) const {
183 return fval_;
184 }
185
186 ROL::Ptr<Vector<Real> > getGradient(void) const {
187 return g_;
188 }
189
191 return nfval_;
192 }
193
195 return ngval_;
196 }
197
198 ROL::Ptr<const Vector<Real> > getLowerMask(void) const {
199 return maskL_;
200 }
201
202 ROL::Ptr<const Vector<Real> > getUpperMask(void) const {
203 return maskU_;
204 }
205
213 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
214 obj_->update(x,flag,iter);
215 }
216
228 Real value( const Vector<Real> &x, Real &tol ) {
229
231 Elementwise::ReductionSum<Real> sum;
232 Elementwise::Multiply<Real> mult;
233
234 // Compute the unpenalized objective value
235 fval_ = obj_->value(x,tol);
236 nfval_++;
237
238 Real fval = fval_;
239 Real linearTerm = 0.0; // Linear damping contribution
240
241 a_->set(x); // a_i = x_i
242 a_->axpy(-1.0,*lo_); // a_i = x_i-l_i
243
244 if( useLinearDamping_ ) {
245
246 c_->set(*maskL_); // c_i = { 1 if l_i > NINF
247 // { 0 otherwise
248 c_->applyBinary(Mask(true),*maskU_); // c_i = { 1 if l_i > NINF and u_i = INF
249 // { 0 otherwise
250 c_->applyBinary(mult,*a_); // c_i = { x_i-l_i if l_i > NINF and u_i = INF
251
252 // Penalizes large positive x_i when only a lower bound exists
253 linearTerm += c_->reduce(sum);
254 }
255
256 a_->applyUnary(mlog); // a_i = mlog(x_i-l_i)
257
258 Real aval = a_->dot(*maskL_);
259
260 b_->set(*up_); // b_i = u_i
261 b_->axpy(-1.0,x); // b_i = u_i-x_i
262
263 if( useLinearDamping_ ) {
264
265 c_->set(*maskU_); // c_i = { 1 if u_i < INF
266 // { 0 otherwise
267 c_->applyBinary(Mask(true),*maskL_); // c_i = { 1 if u_i < INF and l_i = NINF
268 // { 0 otherwise
269 c_->applyBinary(mult,*b_); // c_i = { u_i-x_i if u_i < INF and l_i = NINF
270 // { 0 otherwise
271
272 // Penalizes large negative x_i when only an upper bound exists
273 linearTerm += c_->reduce(sum);
274
275 }
276
277 b_->applyUnary(mlog); // b_i = mlog(u_i-x_i)
278
279 Real bval = b_->dot(*maskU_);
280
281
282 fval -= mu_*(aval+bval);
283 fval += kappaD_*mu_*linearTerm;
284
285 return fval;
286
287 }
288
297 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
298 // Compute gradient of objective function
299 obj_->gradient(*g_,x,tol);
300 ngval_++;
301 g.set(*g_);
302
303 // Add gradient of the log barrier penalty
305
306 a_->set(x); // a = x
307 a_->axpy(-1.0,*lo_); // a = x-l
308
309 a_->applyUnary(mrec); // a_i = 1/(x_i-l_i) for i s.t. x_i > l_i, 0 otherwise
310 a_->applyBinary(Mask(true),*maskL_); // zero elements where l = NINF
311
312 b_->set(*up_); // b = u
313 b_->axpy(-1.0,x); // b = u-x
314 b_->applyUnary(mrec); // b_i = 1/(u_i-x_i) for i s.t. u_i > x_i, 0 otherwise
315 b_->applyBinary(Mask(true),*maskU_); // zero elements where u = INF
316
317 g.axpy(-mu_,*a_);
318 g.axpy(mu_,*b_);
319
320 if( useLinearDamping_ ) {
321
322 a_->set(*maskL_);
323 a_->applyBinary(Mask(true),*maskU_); // a_i = { 1 if l_i > NINF and u_i = INF
324 // { 0 otherwise
325 b_->set(*maskU_);
326 b_->applyBinary(Mask(true),*maskL_); // b_i = { 1 if u_i < INF and l_i = NINF
327 // { 0 otherwise
328 g.axpy(-mu_*kappaD_,*a_);
329 g.axpy( mu_*kappaD_,*b_);
330
331 }
332 }
333
343 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
344
346 Elementwise::Multiply<Real> mult;
347 Elementwise::Power<Real> square(2.0);
348
349 obj_->hessVec(hv,v,x,tol);
350
351 a_->set(x); // a = x
352 a_->axpy(-1.0,*lo_); // a = x-l
353 a_->applyUnary(mrec); // a_i = 1/(x_i-l_i) for i s.t. x_i > l_i, 0 otherwise
354 a_->applyBinary(Mask(true),*maskL_); // zero elements where l = NINF
355 a_->applyUnary(square); // a_i = { (x_i-l_i)^(-2) if l_i > NINF
356 // { 0 if l_i = NINF
357 a_->applyBinary(mult,v);
358
359 b_->set(*up_); // b = u
360 b_->axpy(-1.0,x); // b = u-x
361 b_->applyUnary(mrec); // b_i = 1/(u_i-x_i) for i s.t. u_i > x_i, 0 otherwise
362 b_->applyBinary(Mask(true),*maskU_); // zero elements where u = INF
363 b_->applyUnary(square); // b_i = { (u_i-x_i)^(-2) if u_i < INF
364 // { 0 if u_i = INF
365 b_->applyBinary(mult,v);
366
367 hv.axpy(mu_,*a_);
368 hv.axpy(mu_,*b_);
369 }
370
371 // Return the unpenalized objective
372 const ROL::Ptr<OBJ> getObjective( void ) { return obj_; }
373
374 // Return the bound constraint
375 const ROL::Ptr<BND> getBoundConstraint( void ) { return bnd_; }
376
377}; // class InteriorPointPenalty
378
379}
380
381
382#endif // ROL_INTERIORPOINTPENALTY_H
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Elementwise::ValueSet< Real > ValueSet
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update the interior point penalized objective.
InteriorPointPenalty(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &con, ROL::ParameterList &parlist)
ROL::Ptr< Vector< Real > > getGradient(void) const
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
const ROL::Ptr< OBJ > getObjective(void)
ROL::Ptr< const Vector< Real > > getUpperMask(void) const
const ROL::Ptr< BND > getBoundConstraint(void)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Compute action of Hessian on vector.
ROL::Ptr< const Vector< Real > > getLowerMask(void) const
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Provides the interface to evaluate objective functions.
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .