ROL
ROL_KelleySachsModel.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_KELLEYSACHSMODEL_HPP
11 #define ROL_KELLEYSACHSMODEL_HPP
12 
13 #include "ROL_TrustRegionModel.hpp"
14 #include "ROL_BoundConstraint.hpp"
15 
24 namespace ROL {
25 
26 template<class Real>
27 class KelleySachsModel : public TrustRegionModel<Real> {
28 private:
29  Ptr<Vector<Real>> dual_, prim_, prim2_;
30  Real eps_;
31 
32  class UpperBinding : public Elementwise::BinaryFunction<Real> {
33  public:
34  UpperBinding(Real offset) : offset_(offset) {}
35  Real apply( const Real &x, const Real &y ) const {
36  const Real one(1), tol(1e2*ROL_EPSILON<Real>());
37  return ((y <= -offset_ && x <= tol) ? -one : one);
38  }
39  private:
40  Real offset_;
41  };
42 
43  class LowerBinding : public Elementwise::BinaryFunction<Real> {
44  public:
45  LowerBinding(Real offset) : offset_(offset) {}
46  Real apply( const Real &x, const Real &y ) const {
47  const Real one(1), tol(1e2*ROL_EPSILON<Real>());
48  return ((y >= offset_ && x <= tol) ? -one : one);
49  }
50  private:
51  Real offset_;
52  };
53 
54  class PruneBinding : public Elementwise::BinaryFunction<Real> {
55  public:
56  Real apply( const Real &x, const Real &y ) const {
57  const Real zero(0), one(1);
58  return ((y == one) ? x : zero);
59  }
60  } binding_;
61 
62  class PruneNonbinding : public Elementwise::BinaryFunction<Real> {
63  public:
64  Real apply( const Real &x, const Real &y ) const {
65  const Real zero(0), one(1);
66  return ((y == -one) ? x : zero);
67  }
68  } nonbinding_;
69 
71  const Ptr<const Vector<Real>> gc = TrustRegionModel<Real>::getGradient();
72  const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
73  //const Real one(1);
74  //if (TrustRegionModel<Real>::getBoundConstraint()->isLowerActivated()) {
75  // prim2_->set(*xc);
76  // prim2_->axpy(-one,*TrustRegionModel<Real>::getBoundConstraint()->getLowerBound());
77  // LowerBinding op(eps_);
78  // prim2_->applyBinary(op,*gc);
79  // v.applyBinary(binding_,*prim2_);
80  //}
81  //if (TrustRegionModel<Real>::getBoundConstraint()->isUpperActivated()) {
82  // prim2_->set(*TrustRegionModel<Real>::getBoundConstraint()->getUpperBound());
83  // prim2_->axpy(-one,*xc);
84  // UpperBinding op(eps_);
85  // prim2_->applyBinary(op,*gc);
86  // v.applyBinary(binding_,*prim2_);
87  //}
88  TrustRegionModel<Real>::getBoundConstraint()->pruneActive(v,*gc,*xc,eps_);
89  }
90 
92  const Ptr<const Vector<Real>> gc = TrustRegionModel<Real>::getGradient();
93  const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
94  //const Real one(1);
95  //if (TrustRegionModel<Real>::getBoundConstraint()->isLowerActivated()) {
96  // prim2_->set(*xc);
97  // prim2_->axpy(-one,*TrustRegionModel<Real>::getBoundConstraint()->getLowerBound());
98  // LowerBinding op(eps_);
99  // prim2_->applyBinary(op,*gc);
100  // v.applyBinary(nonbinding_,*prim2_);
101  //}
102  //if (TrustRegionModel<Real>::getBoundConstraint()->isUpperActivated()) {
103  // prim2_->set(*TrustRegionModel<Real>::getBoundConstraint()->getUpperBound());
104  // prim2_->axpy(-one,*xc);
105  // UpperBinding op(eps_);
106  // prim2_->applyBinary(op,*gc);
107  // v.applyBinary(nonbinding_,*prim2_);
108  //}
109  TrustRegionModel<Real>::getBoundConstraint()->pruneInactive(v,*gc,*xc,eps_);
110  }
111 
112 public:
113 
115  const Vector<Real> &x, const Vector<Real> &g,
116  const Ptr<Secant<Real>> &secant = nullPtr,
117  const bool useSecantPrecond = false, const bool useSecantHessVec = false)
118  : TrustRegionModel<Real>::TrustRegionModel(obj,bnd,x,g,secant,useSecantPrecond,useSecantHessVec),
119  eps_(1) {
120  prim_ = x.clone();
121  dual_ = g.clone();
122  prim2_ = x.clone();
123  }
124 
125  void setEpsilon(const Real eps) {
126  eps_ = eps;
127  }
128 
129  /***************************************************************************/
130  /********* BEGIN OBJECTIVE FUNCTION DEFINITIONS **************************/
131  /***************************************************************************/
132  Real value( const Vector<Real> &s, Real &tol ) {
133  hessVec(*dual_,s,s,tol);
134  dual_->scale(static_cast<Real>(0.5));
135  // Remove active components of gradient
138  // Add reduced gradient to reduced hessian in direction s
139  dual_->plus(prim_->dual());
140  return dual_->dot(s.dual());
141  }
142 
143  void gradient( Vector<Real> &g, const Vector<Real> &s, Real &tol ) {
144  // Apply (reduced) hessian to direction s
145  hessVec(g,s,s,tol);
146  // Remove active components of gradient
149  // Add reduced gradient to reduced hessian in direction s
150  g.plus(prim_->dual());
151  }
152 
153  void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
154  // Set vnew to v
155  prim_->set(v);
156  // Remove elements of vnew corresponding to binding set
158  // Apply full Hessian to reduced vector
160  // Remove elements of Hv corresponding to binding set
162  // Set vnew to v
163  prim_->set(v);
164  // Remove Elements of vnew corresponding to complement of binding set
166  dual_->set(prim_->dual());
168  // Fill complement of binding set elements in Hp with v
169  Hv.plus(*dual_);
170  }
171 
172  void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
173  // Set vnew to v
174  dual_->set(v);
175  // Remove elements of vnew corresponding to binding set
177  // Apply full Hessian to reduced vector
179  // Remove elements of Hv corresponding to binding set
181  // Set vnew to v
182  dual_->set(v);
183  // Remove Elements of vnew corresponding to complement of binding set
185  prim_->set(dual_->dual());
187  // Fill complement of binding set elements in Hv with v
188  Hv.plus(*prim_);
189  }
190 
191  void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
192  // Set vnew to v
193  dual_->set(v);
194  // Remove elements of vnew corresponding to binding set
196  // Apply full Hessian to reduced vector
198  // Remove elements of Mv corresponding to binding set
200  // Set vnew to v
201  dual_->set(v);
202  // Remove Elements of vnew corresponding to complement of binding set
204  prim_->set(dual_->dual());
206  // Fill complement of binding set elements in Mv with v
207  Mv.plus(*prim_);
208  }
209  /***************************************************************************/
210  /********* END OBJECTIVE FUNCTION DEFINITIONS ****************************/
211  /***************************************************************************/
212 
213  void dualTransform( Vector<Real> &tv, const Vector<Real> &v ) {
214  // Compute T(v) = P_I(v) where P_I is the projection onto the inactive indices
215  tv.set(v);
217  }
218 
219  void primalTransform( Vector<Real> &tv, const Vector<Real> &v ) {
220  // Compute T(v) = P( x + v ) - x where P is the projection onto the feasible set
221  const Ptr<const Vector<Real>> xc = TrustRegionModel<Real>::getIterate();
222  tv.set(*xc);
223  tv.plus(v);
225  tv.axpy(static_cast<Real>(-1),*xc);
226  }
227 
228 };
229 
230 } // namespace ROL
231 
232 #endif
Provides the interface to evaluate objective functions.
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:192
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Real apply(const Real &x, const Real &y) const
virtual void plus(const Vector &x)=0
Compute , where .
Real apply(const Real &x, const Real &y) const
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
void precond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
void gradient(Vector< Real > &g, const Vector< Real > &s, Real &tol)
Compute gradient.
Provides the interface to evaluate trust-region model functions.
ROL::KelleySachsModel::PruneNonbinding nonbinding_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual const Ptr< const Vector< Real > > getGradient(void) const
ROL::KelleySachsModel::PruneBinding binding_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void pruneNonbindingConstraints(Vector< Real > &v)
void applyPrecond(Vector< Real > &Pv, const Vector< Real > &v, Real &tol)
void applyHessian(Vector< Real > &hv, const Vector< Real > &v, Real &tol)
Real apply(const Real &x, const Real &y) const
void applyInvHessian(Vector< Real > &hv, const Vector< Real > &v, Real &tol)
void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:45
Real apply(const Real &x, const Real &y) const
Ptr< Vector< Real > > dual_
Provides the interface to apply upper and lower bound constraints.
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
Ptr< Vector< Real > > prim2_
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply inverse Hessian approximation to vector.
void setEpsilon(const Real eps)
Provides the interface to evaluate projected trust-region model functions from the Kelley-Sachs bound...
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 pruneBindingConstraints(Vector< Real > &v)
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
Real value(const Vector< Real > &s, Real &tol)
Compute value.
virtual const Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Ptr< Vector< Real > > prim_
void primalTransform(Vector< Real > &tv, const Vector< Real > &v)
virtual const Ptr< const Vector< Real > > getIterate(void) const