ROL
ROL_ObjectiveMMA.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_OBJECTIVEMMA_H
11 #define ROL_OBJECTIVEMMA_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_BoundConstraint.hpp"
15 
25 namespace ROL {
26 
27 template<class Real>
28 class ObjectiveMMA : public Objective<Real> {
29 
30  template <typename T> using ROL::Ptr = ROL::Ptr<T>;
31 
34 
35 private:
36 
37  const ROL::Ptr<OBJ> obj_;
38  const ROL::Ptr<BND> bnd_;
39 
40 
41  ROL::Ptr<V> l_; // Lower bound
42  ROL::Ptr<V> u_; // Upper bound
43 
44  ROL::Ptr<V> p_; // First MMA numerator
45  ROL::Ptr<V> q_; // Second MMA numerator
46 
47  ROL::Ptr<V> d_; // Scratch vector
48 
49  Real fval_; // Original objective value
50 
51  Real tol_;
52 
53 public:
54 
55  ObjectiveMMA( const ROL::Ptr<Objective<Real> > &obj,
56  const ROL::Ptr<BoundConstraint<Real> > &bnd,
57  const Vector<Real> &x,
58  Real tol=std::sqrt(ROL_EPSILON<Real>()) ) :
59  obj_(obj), bnd_(bnd), tol_(tol) {
60 
61  l_ = bnd_->getLowerBound();
62  u_ = bnd_->getUpperBound();
63 
64  p_ = x.clone();
65  q_ = x.clone();
66  d_ = x.clone();
67 
68  }
69 
70  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
71 
72  Elementwise::ThresholdUpper<Real> positive(0.0);
73  Elementwise::Power<Real> square(2.0);
74  Elementwise::Multiply<Real> mult;
75 
76  obj_->update(x,flag,iter);
77 
78  fval_ = obj_->value(x,tol);
79  obj_->gradient(*p_,x,tol);
80  q_->set(*p_);
81 
82  p_->applyUnary(positive);
83  q_->applyUnary(negative);
84 
85  d_->set(x);
86  d_->axpy(-1.0,*l_);
87  d_->applyUnary(square);
88  p_->applyBinary(mult,*d_);
89 
90  d_->set(*u_);
91  d_->axpy(-1.0,x);
92  d_->applyUnary(square);
93  q_->applyBinary(mult,*d_);
94 
95  }
96 
97  /*
98  \f[ F(x) \approx F(x^0) + \sum\limit_{i=1}^n \left( \frac{p_i}{U_i-x_i} + \frac{q_i}{x_i-L_i}\right) \f]
99  */
100  Real value( const Vector<Real> &x, Real &tol ) {
101 
102  Elementwise::ReductionSum<Real> sum;
103  Elementwise::DivideAndInvert<Real> divinv;
104  Real fval = fval_;
105 
106  d_->set(*u_);
107  d_->axpy(-1.0,x);
108  d_->applyBinary(divinv,*p_);
109 
110  fval += d_->reduce(sum);
111 
112  d_->set(x);
113  d_->axpy(-1.0,*l_);
114  d_->applyBinary(divinv,*q_);
115 
116  fval += d_->reduce(sum);
117 
118  return fval;
119 
120  }
121 
122  /*
123  \f[ \frac{F(x)}{\partial x_j} = \frac{p_j}{(U_j-x_j)^2} - \frac{q_j}({x_j-L_j)^2}\ \f]
124  */
125  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
126 
127  Elementwise::DivideAndInvert<Real> divinv;
128  Elementwise::Power<Real> square(2.0);
129 
130  d_->set(*u_);
131  d_->axpy(-1.0,x);
132  d_->applyUnary(square);
133  d_->applyBinary(divinv,*p_);
134 
135  g.set(*d_);
136 
137  d_->set(x);
138  d_->axpy(-1.0,*l_);
139  d_->applyUnary(square);
140  d_->applyBinary(divinv,*q_);
141 
142  g.plus(*d_);
143 
144  }
145 
146  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
147 
148  Elementwise::DivideAndInvert<Real> divinv;
149  Elementwise::Multiply<Real> mult;
150  Elementwise::Power<Real> cube(3.0);
151 
152  d_->set(*u_);
153  d_->axpy(-1.0,x);
154  d_->applyUnary(cube);
155  d_->applyBinary(divinv,*p_);
156  d_->scale(-2.0);
157 
158  hv.set(*d_);
159 
160  d_->set(x);
161  d_->axpy(-1.0,*l_);
162  d_->applyUnary(cube);
163  d_->applyBinary(divinv,*q_);
164  d_->scale(2.0);
165 
166  hv.plus(*d_);
167  hv.applyBinary(mult,v);
168 
169  }
170 
171  void invHessVec( Vector<Real> &h, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
172 
173  Elementwise::DivideAndInvert<Real> divinv;
174  Elementwise::Multiply<Real> mult;
175  Elementwise::Power<Real> cube(3.0);
176 
177  d_->set(*u_);
178  d_->axpy(-1.0,x);
179  d_->applyUnary(cube);
180  d_->applyBinary(divinv,*p_);
181  d_->scale(-2.0);
182 
183  hv.set(*d_);
184 
185  d_->set(x);
186  d_->axpy(-1.0,*l_);
187  d_->applyUnary(cube);
188  d_->applyBinary(divinv,*q_);
189  d_->scale(2.0);
190 
191  hv.plus(*d_);
192  hv.applyBinary(divinv,v);
193 
194  }
195 
196 }; // class ObjectiveMMA
197 
198 } // namespace ROL
199 
200 
201 
202 
203 
204 #endif // ROL_OBJECTIVEMMA_H
205 
Provides the interface to evaluate objective functions.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void plus(const Vector &x)=0
Compute , where .
Provides the interface to to Method of Moving Asymptotes Objective function.
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:214
BoundConstraint< Real > BND
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ObjectiveMMA(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const Vector< Real > &x, Real tol=std::sqrt(ROL_EPSILON< Real >()))
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void invHessVec(Vector< Real > &h, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
const ROL::Ptr< OBJ > obj_
Objective< Real > OBJ
Provides the interface to apply upper and lower bound constraints.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
const ROL::Ptr< BND > bnd_