ROL
ROL_MoreauYosidaPenalty.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_MOREAUYOSIDAPENALTY_H
45 #define ROL_MOREAUYOSIDAPENALTY_H
46 
47 #include "ROL_Objective.hpp"
48 #include "ROL_BoundConstraint.hpp"
49 #include "ROL_Vector.hpp"
50 #include "ROL_Types.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include <iostream>
53 
62 namespace ROL {
63 
64 template <class Real>
65 class MoreauYosidaPenalty : public Objective<Real> {
66 private:
67  const Teuchos::RCP<Objective<Real> > obj_;
68  const Teuchos::RCP<BoundConstraint<Real> > con_;
69 
70  Teuchos::RCP<Vector<Real> > g_;
71  Teuchos::RCP<Vector<Real> > l_;
72  Teuchos::RCP<Vector<Real> > u_;
73  Teuchos::RCP<Vector<Real> > l1_;
74  Teuchos::RCP<Vector<Real> > u1_;
75  Teuchos::RCP<Vector<Real> > dl1_;
76  Teuchos::RCP<Vector<Real> > du1_;
77  Teuchos::RCP<Vector<Real> > xlam_;
78  Teuchos::RCP<Vector<Real> > v_;
79  Teuchos::RCP<Vector<Real> > dv_;
80  Teuchos::RCP<Vector<Real> > dv2_;
81  Teuchos::RCP<Vector<Real> > lam_;
82  Teuchos::RCP<Vector<Real> > tmp_;
83 
84  Real mu_;
85  Real fval_;
87  int nfval_;
88  int ngval_;
89 
90  void computePenalty(const Vector<Real> &x) {
91  if ( con_->isActivated() ) {
92  Real one = 1.0;
93  if ( !isConEvaluated_ ) {
94  xlam_->set(x);
95  xlam_->axpy(one/mu_,*lam_);
96 
97  if ( con_->isFeasible(*xlam_) ) {
98  l1_->zero(); dl1_->zero();
99  u1_->zero(); du1_->zero();
100  }
101  else {
102  // Compute lower penalty component
103  l1_->set(*l_);
104  con_->pruneLowerInactive(*l1_,*xlam_);
105  tmp_->set(*xlam_);
106  con_->pruneLowerInactive(*tmp_,*xlam_);
107  l1_->axpy(-one,*tmp_);
108 
109  // Compute upper penalty component
110  u1_->set(*xlam_);
111  con_->pruneUpperInactive(*u1_,*xlam_);
112  tmp_->set(*u_);
113  con_->pruneUpperInactive(*tmp_,*xlam_);
114  u1_->axpy(-one,*tmp_);
115 
116  // Compute derivative of lower penalty component
117  dl1_->set(l1_->dual());
118  con_->pruneLowerInactive(*dl1_,*xlam_);
119 
120  // Compute derivative of upper penalty component
121  du1_->set(u1_->dual());
122  con_->pruneUpperInactive(*du1_,*xlam_);
123  }
124 
125  isConEvaluated_ = true;
126  }
127  }
128  }
129 
130 public:
132 
133  MoreauYosidaPenalty(const Teuchos::RCP<Objective<Real> > &obj,
134  const Teuchos::RCP<BoundConstraint<Real> > &con,
135  const ROL::Vector<Real> &x, const Real mu = 1.0)
136  : obj_(obj), con_(con), mu_(mu),
137  fval_(0), isConEvaluated_(false), nfval_(0), ngval_(0) {
138 
139  g_ = x.dual().clone();
140  l_ = x.clone();
141  l1_ = x.clone();
142  dl1_ = x.dual().clone();
143  u_ = x.clone();
144  u1_ = x.clone();
145  du1_ = x.dual().clone();
146  xlam_ = x.clone();
147  v_ = x.clone();
148  dv_ = x.dual().clone();
149  dv2_ = x.dual().clone();
150  lam_ = x.clone();
151  tmp_ = x.clone();
152 
153  con_->setVectorToLowerBound(*l_);
154  con_->setVectorToUpperBound(*u_);
155 
156  lam_->zero();
157  //lam_->set(*u_);
158  //lam_->plus(*l_);
159  //lam_->scale(0.5);
160  }
161 
162  void updateMultipliers(Real mu, const ROL::Vector<Real> &x) {
163  if ( con_->isActivated() ) {
164  Real one = 1.0;
165  computePenalty(x);
166 
167  lam_->set(*u1_);
168  lam_->axpy(-one,*l1_);
169  lam_->scale(mu_);
170 
171  mu_ = mu;
172  }
173 
174  nfval_ = 0; ngval_ = 0;
175  isConEvaluated_ = false;
176  }
177 
178  Real getObjectiveValue(void) const {
179  return fval_;
180  }
181 
182  Teuchos::RCP<Vector<Real> > getGradient(void) const {
183  return g_;
184  }
185 
187  return nfval_;
188  }
189 
191  return ngval_;
192  }
193 
201  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
202  obj_->update(x,flag,iter);
203  con_->update(x,flag,iter);
204  isConEvaluated_ = false;
205  }
206 
213  Real value( const Vector<Real> &x, Real &tol ) {
214  Real half = 0.5;
215  // Compute objective function value
216  fval_ = obj_->value(x,tol);
217  nfval_++;
218  // Add value of the Moreau-Yosida penalty
219  Real fval = fval_;
220  if ( con_->isActivated() ) {
221  computePenalty(x);
222  fval += half*mu_*(l1_->dot(*l1_) + u1_->dot(*u1_));
223  }
224  return fval;
225  }
226 
234  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
235  // Compute gradient of objective function
236  obj_->gradient(*g_,x,tol);
237  ngval_++;
238  g.set(*g_);
239  // Add gradient of the Moreau-Yosida penalty
240  if ( con_->isActivated() ) {
241  computePenalty(x);
242  g.axpy(-mu_,*dl1_);
243  g.axpy(mu_,*du1_);
244  }
245  }
246 
255  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
256  // Apply objective Hessian to a vector
257  obj_->hessVec(hv,v,x,tol);
258  // Add Hessian of the Moreau-Yosida penalty
259  if ( con_->isActivated() ) {
260  Real one = 1.0;
261  computePenalty(x);
262 
263  v_->set(v);
264  con_->pruneLowerActive(*v_,*xlam_);
265  v_->scale(-one);
266  v_->plus(v);
267  dv_->set(v_->dual());
268  dv2_->set(*dv_);
269  con_->pruneLowerActive(*dv_,*xlam_);
270  dv_->scale(-one);
271  dv_->plus(*dv2_);
272  hv.axpy(mu_,*dv_);
273 
274  v_->set(v);
275  con_->pruneUpperActive(*v_,*xlam_);
276  v_->scale(-one);
277  v_->plus(v);
278  dv_->set(v_->dual());
279  dv2_->set(*dv_);
280  con_->pruneUpperActive(*dv_,*xlam_);
281  dv_->scale(-one);
282  dv_->plus(*dv2_);
283  hv.axpy(mu_,*dv_);
284  }
285  }
286 
287 // Definitions for parametrized (stochastic) objective functions
288 public:
289  void setParameter(const std::vector<Real> &param) {
291  obj_->setParameter(param);
292  }
293 }; // class MoreauYosidaPenalty
294 
295 } // namespace ROL
296 
297 #endif
Provides the interface to evaluate objective functions.
const Teuchos::RCP< Objective< Real > > obj_
Teuchos::RCP< Vector< Real > > g_
void setParameter(const std::vector< Real > &param)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
Teuchos::RCP< Vector< Real > > dv_
Teuchos::RCP< Vector< Real > > xlam_
Contains definitions of custom data types in ROL.
Teuchos::RCP< Vector< Real > > v_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void updateMultipliers(Real mu, const ROL::Vector< Real > &x)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< Vector< Real > > dv2_
Teuchos::RCP< Vector< Real > > tmp_
void computePenalty(const Vector< Real > &x)
Teuchos::RCP< Vector< Real > > u_
Teuchos::RCP< Vector< Real > > l_
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:213
Teuchos::RCP< Vector< Real > > du1_
Teuchos::RCP< Vector< Real > > getGradient(void) const
MoreauYosidaPenalty(const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< BoundConstraint< Real > > &con, const ROL::Vector< Real > &x, const Real mu=1.0)
Provides the interface to evaluate the Moreau-Yosida penalty function.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Provides the interface to apply upper and lower bound constraints.
const Teuchos::RCP< BoundConstraint< Real > > con_
Teuchos::RCP< Vector< Real > > lam_
virtual void setParameter(const std::vector< Real > &param)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update Moreau-Yosida penalty function.
Teuchos::RCP< Vector< Real > > u1_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< Vector< Real > > dl1_
Teuchos::RCP< Vector< Real > > l1_