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