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 
47 #include "ROL_TrustRegionModel.hpp"
48 #include "ROL_BoundConstraint.hpp"
49 
58 namespace ROL {
59 
60 template<class Real>
61 class KelleySachsModel : public TrustRegionModel<Real> {
62 private:
63  Teuchos::RCP<BoundConstraint<Real> > bnd_;
64  Teuchos::RCP<Secant<Real> > secant_;
65  Teuchos::RCP<Vector<Real> > dual_, prim_;
66 
67  const bool useSecantPrecond_;
68  const bool useSecantHessVec_;
69 
70  Real eps_;
71 
72 public:
73 
75  const Vector<Real> &x, const Vector<Real> &g, const Real eps)
76  : TrustRegionModel<Real>::TrustRegionModel(obj,x,g,false),
77  secant_(Teuchos::null), useSecantPrecond_(false), useSecantHessVec_(false), eps_(eps) {
78  bnd_ = Teuchos::rcpFromRef(bnd);
79  prim_ = x.clone();
80  dual_ = g.clone();
81  }
82 
84  const Vector<Real> &x, const Vector<Real> &g, const Real eps,
85  const Teuchos::RCP<Secant<Real> > &secant,
86  const bool useSecantPrecond, const bool useSecantHessVec)
87  : TrustRegionModel<Real>::TrustRegionModel(obj,x,g,false),
88  secant_(secant), useSecantPrecond_(useSecantPrecond), useSecantHessVec_(useSecantHessVec), eps_(eps) {
89  bnd_ = Teuchos::rcpFromRef(bnd);
90  prim_ = x.clone();
91  dual_ = g.clone();
92  }
93 
94  Real value( const Vector<Real> &s, Real &tol ) {
95  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
96  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
97  hessVec(*dual_,s,s,tol);
98  dual_->scale(static_cast<Real>(0.5));
99  // Remove active components of gradient
100  prim_->set(gc->dual());
101  bnd_->pruneActive(*prim_,*gc,*xc,eps_);
102  // Add reduced gradient to reduced hessian in direction s
103  dual_->plus(prim_->dual());
104  return dual_->dot(s.dual());
105  }
106 
107  void gradient( Vector<Real> &g, const Vector<Real> &s, Real &tol ) {
108  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
109  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
110  // Apply (reduced) hessian to direction s
111  hessVec(g,s,s,tol);
112  // Remove active components of gradient
114  bnd_->pruneActive(*prim_,*gc,*xc,eps_);
115  // Add reduced gradient to reduced hessian in direction s
116  g.plus(prim_->dual());
117  }
118 
119  void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
120  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
121  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
122  // Set vnew to v
123  prim_->set(v);
124  // Remove elements of vnew corresponding to binding set
125  bnd_->pruneActive(*prim_,*gc,*xc,eps_);
126  // Apply full Hessian to reduced vector
127  if ( useSecantHessVec_ ) {
128  secant_->applyB(Hv,*prim_);
129  }
130  else {
131  TrustRegionModel<Real>::getObjective()->hessVec(Hv,*prim_,*xc,tol);
132  }
133  // Remove elements of Hv corresponding to binding set
134  bnd_->pruneActive(Hv,*gc,*xc,eps_);
135  // Set vnew to v
136  prim_->set(v);
137  // Remove Elements of vnew corresponding to complement of binding set
138  bnd_->pruneInactive(*prim_,*gc,*xc,eps_);
139  dual_->set(prim_->dual());
140  bnd_->pruneInactive(*dual_,*gc,*xc,eps_);
141  // Fill complement of binding set elements in Hp with v
142  Hv.plus(*dual_);
143  }
144 
145  void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
146  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
147  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
148  // Set vnew to v
149  dual_->set(v);
150  // Remove elements of vnew corresponding to binding set
151  bnd_->pruneActive(*dual_,*gc,*xc,eps_);
152  // Apply full Hessian to reduced vector
153  if ( useSecantHessVec_ ) {
154  secant_->applyH(Hv,*dual_);
155  }
156  else {
157  TrustRegionModel<Real>::getObjective()->invHessVec(Hv,*dual_,*xc,tol);
158  }
159  // Remove elements of Hv corresponding to binding set
160  bnd_->pruneActive(Hv,*gc,*xc,eps_);
161  // Set vnew to v
162  dual_->set(v);
163  // Remove Elements of vnew corresponding to complement of binding set
164  bnd_->pruneInactive(*dual_,*gc,*xc,eps_);
165  prim_->set(dual_->dual());
166  bnd_->pruneInactive(*prim_,*gc,*xc,eps_);
167  // Fill complement of binding set elements in Hv with v
168  Hv.plus(*prim_);
169  }
170 
171  void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
172  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
173  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
174  // Set vnew to v
175  dual_->set(v);
176  // Remove elements of vnew corresponding to binding set
177  bnd_->pruneActive(*dual_,*gc,*xc,eps_);
178  // Apply full Hessian to reduced vector
179  if ( useSecantPrecond_ ) {
180  secant_->applyH(Mv,*dual_);
181  }
182  else {
183  TrustRegionModel<Real>::getObjective()->precond(Mv,*dual_,*xc,tol);
184  }
185  // Remove elements of Mv corresponding to binding set
186  bnd_->pruneActive(Mv,*gc,*xc,eps_);
187  // Set vnew to v
188  dual_->set(v);
189  // Remove Elements of vnew corresponding to complement of binding set
190  bnd_->pruneInactive(*dual_,*gc,*xc,eps_);
191  prim_->set(dual_->dual());
192  bnd_->pruneInactive(*prim_,*gc,*xc,eps_);
193  // Fill complement of binding set elements in Mv with v
194  Mv.plus(*prim_);
195  }
196 
197  void dualTransform( Vector<Real> &tv, const Vector<Real> &v ) {
198  // Compute T(v) = P_I(v) where P_I is the projection onto the inactive indices
199  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
200  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
201  tv.set(v);
202  bnd_->pruneActive(tv,*gc,*xc,eps_);
203  }
204 
205  void primalTransform( Vector<Real> &tv, const Vector<Real> &v ) {
206  // Compute T(v) = P( x + v ) - x where P is the projection onto the feasible set
207  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
208  tv.set(*xc);
209  tv.plus(v);
210  bnd_->project(tv);
211  tv.axpy(static_cast<Real>(-1),*xc);
212  }
213 
214  const Teuchos::RCP<BoundConstraint<Real> > getBoundConstraint(void) const {
215  return bnd_;
216  }
217 
218 };
219 
220 } // namespace ROL
221 
222 #endif
Provides the interface to evaluate objective functions.
KelleySachsModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g, const Real eps, const Teuchos::RCP< Secant< Real > > &secant, const bool useSecantPrecond, const bool useSecantHessVec)
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
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.
const Teuchos::RCP< BoundConstraint< Real > > getBoundConstraint(void) const
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Provides the interface to evaluate trust-region model functions.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
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 > > prim_
KelleySachsModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g, const Real eps)
Teuchos::RCP< BoundConstraint< Real > > bnd_
void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
virtual const Teuchos::RCP< const Vector< Real > > getGradient(void) const
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:70
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.
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply inverse Hessian approximation to vector.
Provides the interface to evaluate projected trust-region model functions from the Kelley-Sachs bound...
Teuchos::RCP< Secant< Real > > secant_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
Real value(const Vector< Real > &s, Real &tol)
Compute value.
Teuchos::RCP< Vector< Real > > dual_
void primalTransform(Vector< Real > &tv, const Vector< Real > &v)
virtual const Teuchos::RCP< Objective< Real > > getObjective(void) const
virtual const Teuchos::RCP< const Vector< Real > > getIterate(void) const