ROL
ROL_AugmentedLagrangian.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_AUGMENTEDLAGRANGIAN_H
45 #define ROL_AUGMENTEDLAGRANGIAN_H
46 
47 #include "ROL_Objective.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 AugmentedLagrangian : public Objective<Real> {
66 private:
67  Teuchos::RCP<Objective<Real> > obj_;
68  Teuchos::RCP<EqualityConstraint<Real> > con_;
69  Teuchos::RCP<Vector<Real> > lam_;
70  Teuchos::RCP<Vector<Real> > dlam_;
71  Teuchos::RCP<Vector<Real> > x_;
72  Teuchos::RCP<Vector<Real> > g_;
73  Teuchos::RCP<Vector<Real> > c_;
74  Teuchos::RCP<Vector<Real> > dc1_;
75  Teuchos::RCP<Vector<Real> > dc2_;
76 
77  Real fval_;
78  int ncval_;
79  int nfval_;
80  int ngval_;
81 
85 
89 
90 public:
92 
94  const Vector<Real> &x, const Vector<Real> &c,
95  const Vector<Real> &l, const Real mu,
96  Teuchos::ParameterList &parlist)
97  : fval_(0.0), ncval_(0), nfval_(0), ngval_(0), penaltyParameter_(mu),
98  isConstraintComputed_(false), isValueComputed_(false), isGradientComputed_(false) {
99  obj_ = Teuchos::rcp(&obj, false);
100  con_ = Teuchos::rcp(&con, false);
101 
102  x_ = x.clone();
103  g_ = x.dual().clone();
104  dc1_ = x.dual().clone();
105  dc2_ = c.clone();
106  c_ = c.clone();
107  lam_ = l.clone();
108  dlam_ = l.clone();
109  lam_->set(l);
110 
111  Teuchos::ParameterList& sublist = parlist.sublist("Step").sublist("Augmented Lagrangian");
112  scaleLagrangian_ = sublist.get("Use Scaled Augmented Lagrangian", false);
113  HessianLevel_ = sublist.get("Level of Hessian Approximation", 0);
114  }
115 
117  Real tol = std::sqrt(ROL_EPSILON);
118  if ( !isValueComputed_ ) {
119  // Evaluate objective function value
120  fval_ = obj_->value(x,tol);
121  nfval_++;
122  isValueComputed_ = true;
123  }
124  return fval_;
125  }
126 
128  Real tol = std::sqrt(ROL_EPSILON);
129  if ( !isGradientComputed_ ) {
130  // Compute objective function gradient
131  obj_->gradient(*g_,x,tol);
132  ngval_++;
133  isGradientComputed_ = true;
134  }
135  g.set(*g_);
136  }
137 
139  Real tol = std::sqrt(ROL_EPSILON);
140  if ( !isConstraintComputed_ ) {
141  // Evaluate constraint
142  con_->value(*c_,x,tol);
143  ncval_++;
144  isConstraintComputed_ = true;
145  }
146  c.set(*c_);
147  }
148 
150  return ncval_;
151  }
152 
154  return nfval_;
155  }
156 
158  return ngval_;
159  }
160 
161  void reset(const Vector<Real> &lam, const Real penaltyParameter) {
162  ncval_ = 0.; nfval_ = 0.; ngval_ = 0.;
163  lam_->set(lam);
164  penaltyParameter_ = penaltyParameter;
165  }
166 
174  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
175  obj_->update(x,flag,iter);
176  con_->update(x,flag,iter);
177  if ( flag ) {
178  isConstraintComputed_ = false;
179  isValueComputed_ = false;
180  isGradientComputed_ = false;
181  }
182  }
183 
190  Real value( const Vector<Real> &x, Real &tol ) {
191  if ( !isConstraintComputed_ ) {
192  // Evaluate constraint
193  con_->value(*c_,x,tol);
194  ncval_++;
195  isConstraintComputed_ = true;
196  }
197  if ( !isValueComputed_ ) {
198  // Compute objective function value
199  fval_ = obj_->value(x,tol);
200  nfval_++;
201  isValueComputed_ = true;
202  }
203  // Apply Lagrange multiplier to constraint
204  Real cval = lam_->dot(c_->dual());
205  // Compute penalty term
206  Real pval = c_->dot(*c_);
207  // Compute Augmented Lagrangian value
208  Real val = 0.0;
209  if (scaleLagrangian_) {
210  val = (fval_ + cval)/penaltyParameter_ + 0.5*pval;
211  }
212  else {
213  val = fval_ + cval + 0.5*penaltyParameter_*pval;
214  }
215  return val;
216  }
217 
225  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
226  if ( !isConstraintComputed_ ) {
227  // Evaluate constraint
228  con_->value(*c_,x,tol);
229  ncval_++;
230  isConstraintComputed_ = true;
231  }
232  if ( !isGradientComputed_ ) {
233  // Compute objective function gradient
234  obj_->gradient(*g_,x,tol);
235  ngval_++;
236  isGradientComputed_ = true;
237  }
238  g.set(*g_);
239  // Compute gradient of Augmented Lagrangian
240  dlam_->set(c_->dual());
241  if ( scaleLagrangian_ ) {
242  g.scale(1./penaltyParameter_);
243  dlam_->axpy(1./penaltyParameter_,*lam_);
244  }
245  else {
246  dlam_->scale(penaltyParameter_);
247  dlam_->plus(*lam_);
248  }
249  con_->applyAdjointJacobian(*dc1_,*dlam_,x,tol);
250  g.plus(*dc1_);
251  }
252 
261  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
262  // Apply objective Hessian to a vector
263  obj_->hessVec(hv,v,x,tol);
264  if (HessianLevel_ < 2) {
265  con_->applyJacobian(*dc2_,v,x,tol);
266  con_->applyAdjointJacobian(*dc1_,dc2_->dual(),x,tol);
267  if (scaleLagrangian_) {
268  hv.scale(1./penaltyParameter_);
269  hv.plus(*dc1_);
270  }
271  else {
272  hv.axpy(penaltyParameter_,*dc1_);
273  }
274 
275  if (HessianLevel_ == 0) {
276  if ( !isConstraintComputed_ ) {
277  // Evaluate constraint
278  con_->value(*c_,x,tol);
279  ncval_++;
280  isConstraintComputed_ = true;
281  }
282  // Apply Augmented Lagrangian Hessian to a vector
283  dlam_->set(c_->dual());
284  if ( scaleLagrangian_ ) {
285  dlam_->axpy(1./penaltyParameter_,*lam_);
286  }
287  else {
288  dlam_->scale(penaltyParameter_);
289  dlam_->plus(*lam_);
290  }
291  con_->applyAdjointHessian(*dc1_,*dlam_,v,x,tol);
292  hv.plus(*dc1_);
293  }
294  }
295  }
296 
297 }; // class AugmentedLagrangian
298 
299 } // namespace ROL
300 
301 #endif
Provides the interface to evaluate objective functions.
Provides the interface to evaluate the augmented Lagrangian.
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
virtual void scale(const Real alpha)=0
Compute where .
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update augmented Lagrangian function.
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
Teuchos::RCP< Vector< Real > > g_
Contains definitions of custom data types in ROL.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Teuchos::RCP< Objective< Real > > obj_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< Vector< Real > > c_
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Real getObjectiveValue(const Vector< Real > &x)
Teuchos::RCP< Vector< Real > > dc2_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Defines the equality constraint operator interface.
Teuchos::RCP< EqualityConstraint< Real > > con_
void reset(const Vector< Real > &lam, const Real penaltyParameter)
void getConstraintVec(Vector< Real > &c, const Vector< Real > &x)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Teuchos::RCP< Vector< Real > > dc1_
void getObjectiveGradient(Vector< Real > &g, const Vector< Real > &x)
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
Teuchos::RCP< Vector< Real > > lam_
Teuchos::RCP< Vector< Real > > x_
AugmentedLagrangian(Objective< Real > &obj, EqualityConstraint< Real > &con, const Vector< Real > &x, const Vector< Real > &c, const Vector< Real > &l, const Real mu, Teuchos::ParameterList &parlist)
Teuchos::RCP< Vector< Real > > dlam_
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118