ROL
ROL_KelleySachsModel.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_KELLEYSACHSMODEL_HPP
45#define ROL_KELLEYSACHSMODEL_HPP
46
49
58namespace ROL {
59
60template<class Real>
61class KelleySachsModel : public TrustRegionModel<Real> {
62private:
63 Ptr<Vector<Real>> dual_, prim_, prim2_;
64 Real eps_;
65
66 class UpperBinding : public Elementwise::BinaryFunction<Real> {
67 public:
68 UpperBinding(Real offset) : offset_(offset) {}
69 Real apply( const Real &x, const Real &y ) const {
70 const Real one(1), tol(1e2*ROL_EPSILON<Real>());
71 return ((y <= -offset_ && x <= tol) ? -one : one);
72 }
73 private:
74 Real offset_;
75 };
76
77 class LowerBinding : public Elementwise::BinaryFunction<Real> {
78 public:
79 LowerBinding(Real offset) : offset_(offset) {}
80 Real apply( const Real &x, const Real &y ) const {
81 const Real one(1), tol(1e2*ROL_EPSILON<Real>());
82 return ((y >= offset_ && x <= tol) ? -one : one);
83 }
84 private:
85 Real offset_;
86 };
87
88 class PruneBinding : public Elementwise::BinaryFunction<Real> {
89 public:
90 Real apply( const Real &x, const Real &y ) const {
91 const Real zero(0), one(1);
92 return ((y == one) ? x : zero);
93 }
95
96 class PruneNonbinding : public Elementwise::BinaryFunction<Real> {
97 public:
98 Real apply( const Real &x, const Real &y ) const {
99 const Real zero(0), one(1);
100 return ((y == -one) ? x : zero);
101 }
103
105 const Ptr<const Vector<Real>> gc = TrustRegionModel<Real>::getGradient();
106 const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
107 //const Real one(1);
108 //if (TrustRegionModel<Real>::getBoundConstraint()->isLowerActivated()) {
109 // prim2_->set(*xc);
110 // prim2_->axpy(-one,*TrustRegionModel<Real>::getBoundConstraint()->getLowerBound());
111 // LowerBinding op(eps_);
112 // prim2_->applyBinary(op,*gc);
113 // v.applyBinary(binding_,*prim2_);
114 //}
115 //if (TrustRegionModel<Real>::getBoundConstraint()->isUpperActivated()) {
116 // prim2_->set(*TrustRegionModel<Real>::getBoundConstraint()->getUpperBound());
117 // prim2_->axpy(-one,*xc);
118 // UpperBinding op(eps_);
119 // prim2_->applyBinary(op,*gc);
120 // v.applyBinary(binding_,*prim2_);
121 //}
122 TrustRegionModel<Real>::getBoundConstraint()->pruneActive(v,*gc,*xc,eps_);
123 }
124
126 const Ptr<const Vector<Real>> gc = TrustRegionModel<Real>::getGradient();
127 const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
128 //const Real one(1);
129 //if (TrustRegionModel<Real>::getBoundConstraint()->isLowerActivated()) {
130 // prim2_->set(*xc);
131 // prim2_->axpy(-one,*TrustRegionModel<Real>::getBoundConstraint()->getLowerBound());
132 // LowerBinding op(eps_);
133 // prim2_->applyBinary(op,*gc);
134 // v.applyBinary(nonbinding_,*prim2_);
135 //}
136 //if (TrustRegionModel<Real>::getBoundConstraint()->isUpperActivated()) {
137 // prim2_->set(*TrustRegionModel<Real>::getBoundConstraint()->getUpperBound());
138 // prim2_->axpy(-one,*xc);
139 // UpperBinding op(eps_);
140 // prim2_->applyBinary(op,*gc);
141 // v.applyBinary(nonbinding_,*prim2_);
142 //}
143 TrustRegionModel<Real>::getBoundConstraint()->pruneInactive(v,*gc,*xc,eps_);
144 }
145
146public:
147
149 const Vector<Real> &x, const Vector<Real> &g,
150 const Ptr<Secant<Real>> &secant = nullPtr,
151 const bool useSecantPrecond = false, const bool useSecantHessVec = false)
152 : TrustRegionModel<Real>::TrustRegionModel(obj,bnd,x,g,secant,useSecantPrecond,useSecantHessVec),
153 eps_(1) {
154 prim_ = x.clone();
155 dual_ = g.clone();
156 prim2_ = x.clone();
157 }
158
159 void setEpsilon(const Real eps) {
160 eps_ = eps;
161 }
162
163 /***************************************************************************/
164 /********* BEGIN OBJECTIVE FUNCTION DEFINITIONS **************************/
165 /***************************************************************************/
166 Real value( const Vector<Real> &s, Real &tol ) {
167 hessVec(*dual_,s,s,tol);
168 dual_->scale(static_cast<Real>(0.5));
169 // Remove active components of gradient
172 // Add reduced gradient to reduced hessian in direction s
173 dual_->plus(prim_->dual());
174 return dual_->dot(s.dual());
175 }
176
177 void gradient( Vector<Real> &g, const Vector<Real> &s, Real &tol ) {
178 // Apply (reduced) hessian to direction s
179 hessVec(g,s,s,tol);
180 // Remove active components of gradient
183 // Add reduced gradient to reduced hessian in direction s
184 g.plus(prim_->dual());
185 }
186
187 void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
188 // Set vnew to v
189 prim_->set(v);
190 // Remove elements of vnew corresponding to binding set
192 // Apply full Hessian to reduced vector
194 // Remove elements of Hv corresponding to binding set
196 // Set vnew to v
197 prim_->set(v);
198 // Remove Elements of vnew corresponding to complement of binding set
200 dual_->set(prim_->dual());
202 // Fill complement of binding set elements in Hp with v
203 Hv.plus(*dual_);
204 }
205
206 void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
207 // Set vnew to v
208 dual_->set(v);
209 // Remove elements of vnew corresponding to binding set
211 // Apply full Hessian to reduced vector
213 // Remove elements of Hv corresponding to binding set
215 // Set vnew to v
216 dual_->set(v);
217 // Remove Elements of vnew corresponding to complement of binding set
219 prim_->set(dual_->dual());
221 // Fill complement of binding set elements in Hv with v
222 Hv.plus(*prim_);
223 }
224
225 void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
226 // Set vnew to v
227 dual_->set(v);
228 // Remove elements of vnew corresponding to binding set
230 // Apply full Hessian to reduced vector
232 // Remove elements of Mv corresponding to binding set
234 // Set vnew to v
235 dual_->set(v);
236 // Remove Elements of vnew corresponding to complement of binding set
238 prim_->set(dual_->dual());
240 // Fill complement of binding set elements in Mv with v
241 Mv.plus(*prim_);
242 }
243 /***************************************************************************/
244 /********* END OBJECTIVE FUNCTION DEFINITIONS ****************************/
245 /***************************************************************************/
246
248 // Compute T(v) = P_I(v) where P_I is the projection onto the inactive indices
249 tv.set(v);
251 }
252
254 // Compute T(v) = P( x + v ) - x where P is the projection onto the feasible set
255 const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
256 tv.set(*xc);
257 tv.plus(v);
259 tv.axpy(static_cast<Real>(-1),*xc);
260 }
261
262};
263
264} // namespace ROL
265
266#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Real apply(const Real &x, const Real &y) const
Provides the interface to evaluate projected trust-region model functions from the Kelley-Sachs bound...
void precond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
ROL::KelleySachsModel::PruneBinding binding_
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply inverse Hessian approximation to vector.
void primalTransform(Vector< Real > &tv, const Vector< Real > &v)
KelleySachsModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g, const Ptr< Secant< Real > > &secant=nullPtr, const bool useSecantPrecond=false, const bool useSecantHessVec=false)
void gradient(Vector< Real > &g, const Vector< Real > &s, Real &tol)
Compute gradient.
Ptr< Vector< Real > > prim_
void pruneBindingConstraints(Vector< Real > &v)
ROL::KelleySachsModel::PruneNonbinding nonbinding_
void pruneNonbindingConstraints(Vector< Real > &v)
Real value(const Vector< Real > &s, Real &tol)
Compute value.
Ptr< Vector< Real > > dual_
void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
void setEpsilon(const Real eps)
Ptr< Vector< Real > > prim2_
Provides the interface to evaluate objective functions.
Provides interface for and implements limited-memory secant operators.
Provides the interface to evaluate trust-region model functions.
void applyInvHessian(Vector< Real > &hv, const Vector< Real > &v, Real &tol)
virtual const Ptr< const Vector< Real > > getGradient(void) const
void applyHessian(Vector< Real > &hv, const Vector< Real > &v, Real &tol)
void applyPrecond(Vector< Real > &Pv, const Vector< Real > &v, Real &tol)
virtual const Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
virtual const Ptr< const Vector< Real > > getIterate(void) const
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set 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 ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .