ROL
ROL_QuadraticPenalty.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_QUADRATICPENALTY_H
45#define ROL_QUADRATICPENALTY_H
46
47#include "ROL_Objective.hpp"
48#include "ROL_Constraint.hpp"
49#include "ROL_Vector.hpp"
50#include "ROL_Types.hpp"
51#include "ROL_Ptr.hpp"
52#include <iostream>
53
81namespace ROL {
82
83template <class Real>
84class QuadraticPenalty : public Objective<Real> {
85private:
86 // Required for quadratic penalty definition
87 const ROL::Ptr<Constraint<Real> > con_;
88 ROL::Ptr<Vector<Real> > multiplier_;
90
91 // Auxiliary storage
92 ROL::Ptr<Vector<Real> > primalMultiplierVector_;
93 ROL::Ptr<Vector<Real> > dualOptVector_;
94 ROL::Ptr<Vector<Real> > primalConVector_;
95
96 // Constraint evaluations
97 ROL::Ptr<Vector<Real> > conValue_;
98 Real cscale_;
99
100 // Evaluation counters
102
103 // User defined options
104 const bool useScaling_;
105 const int HessianApprox_;
106
107 // Flags to recompute quantities
109
110 void evaluateConstraint(const Vector<Real> &x, Real &tol) {
111 if ( !isConstraintComputed_ ) {
112 // Evaluate constraint
113 con_->value(*conValue_,x,tol); ncval_++;
115 }
116 }
117
118public:
119 QuadraticPenalty(const ROL::Ptr<Constraint<Real> > &con,
120 const Vector<Real> &multiplier,
121 const Real penaltyParameter,
122 const Vector<Real> &optVec,
123 const Vector<Real> &conVec,
124 const bool useScaling = false,
125 const int HessianApprox = 0 )
126 : con_(con), penaltyParameter_(penaltyParameter), cscale_(1), ncval_(0),
127 useScaling_(useScaling), HessianApprox_(HessianApprox), isConstraintComputed_(false) {
128
129 dualOptVector_ = optVec.dual().clone();
130 primalConVector_ = conVec.clone();
131 conValue_ = conVec.clone();
132 multiplier_ = multiplier.clone();
133 primalMultiplierVector_ = multiplier.clone();
134 }
135
136 void setScaling(const Real cscale = 1) {
137 cscale_ = cscale;
138 }
139
140 virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
141 con_->update(x,flag,iter);
142 isConstraintComputed_ = ((flag || (!flag && iter < 0)) ? false : isConstraintComputed_);
143 }
144
145 virtual Real value( const Vector<Real> &x, Real &tol ) {
146 // Evaluate constraint
147 evaluateConstraint(x,tol);
148 // Apply Lagrange multiplier to constraint
149 Real cval = cscale_*multiplier_->dot(conValue_->dual());
150 // Compute penalty term
151 Real pval = cscale_*cscale_*conValue_->dot(*conValue_);
152 // Compute quadratic penalty value
153 const Real half(0.5);
154 Real val(0);
155 if (useScaling_) {
156 val = cval/penaltyParameter_ + half*pval;
157 }
158 else {
159 val = cval + half*penaltyParameter_*pval;
160 }
161 return val;
162 }
163
164 virtual void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
165 // Evaluate constraint
166 evaluateConstraint(x,tol);
167 // Compute gradient of Augmented Lagrangian
169 if ( useScaling_ ) {
172 }
173 else {
176 }
177 con_->applyAdjointJacobian(g,*primalMultiplierVector_,x,tol);
178 }
179
180 virtual void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
181 // Apply objective Hessian to a vector
182 if (HessianApprox_ < 3) {
183 con_->update(x);
184 con_->applyJacobian(*primalConVector_,v,x,tol);
185 con_->applyAdjointJacobian(hv,primalConVector_->dual(),x,tol);
186 if (!useScaling_) {
188 }
189 else {
191 }
192
193 if (HessianApprox_ == 1) {
194 // Apply Augmented Lagrangian Hessian to a vector
196 if ( useScaling_ ) {
198 }
199 else {
201 }
202 con_->applyAdjointHessian(*dualOptVector_,*primalMultiplierVector_,v,x,tol);
203 hv.plus(*dualOptVector_);
204 }
205
206 if (HessianApprox_ == 0) {
207 // Evaluate constraint
208 evaluateConstraint(x,tol);
209 // Apply Augmented Lagrangian Hessian to a vector
211 if ( useScaling_ ) {
214 }
215 else {
218 }
219 con_->applyAdjointHessian(*dualOptVector_,*primalMultiplierVector_,v,x,tol);
220 hv.plus(*dualOptVector_);
221 }
222 }
223 else {
224 hv.zero();
225 }
226 }
227
228 // Return constraint value
229 virtual void getConstraintVec(Vector<Real> &c, const Vector<Real> &x) {
230 Real tol = std::sqrt(ROL_EPSILON<Real>());
231 // Evaluate constraint
232 evaluateConstraint(x,tol);
233 c.set(*conValue_);
234 }
235
236 // Return total number of constraint evaluations
237 virtual int getNumberConstraintEvaluations(void) const {
238 return ncval_;
239 }
240
241 // Reset with upated penalty parameter
242 virtual void reset(const Vector<Real> &multiplier, const Real penaltyParameter) {
243 ncval_ = 0;
244 multiplier_->set(multiplier);
245 penaltyParameter_ = penaltyParameter;
246 }
247}; // class AugmentedLagrangian
248
249} // namespace ROL
250
251#endif
Contains definitions of custom data types in ROL.
Defines the general constraint operator interface.
Provides the interface to evaluate objective functions.
Provides the interface to evaluate the quadratic constraint penalty.
void setScaling(const Real cscale=1)
virtual Real value(const Vector< Real > &x, Real &tol)
Compute value.
ROL::Ptr< Vector< Real > > multiplier_
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
virtual void reset(const Vector< Real > &multiplier, const Real penaltyParameter)
QuadraticPenalty(const ROL::Ptr< Constraint< Real > > &con, const Vector< Real > &multiplier, const Real penaltyParameter, const Vector< Real > &optVec, const Vector< Real > &conVec, const bool useScaling=false, const int HessianApprox=0)
const ROL::Ptr< Constraint< Real > > con_
virtual int getNumberConstraintEvaluations(void) const
void evaluateConstraint(const Vector< Real > &x, Real &tol)
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
ROL::Ptr< Vector< Real > > primalMultiplierVector_
ROL::Ptr< Vector< Real > > primalConVector_
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
ROL::Ptr< Vector< Real > > dualOptVector_
virtual void getConstraintVec(Vector< Real > &c, const Vector< Real > &x)
ROL::Ptr< Vector< Real > > conValue_
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.