ROL
ROL_Secant.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_SECANT_H
45 #define ROL_SECANT_H
46 
51 #include "Teuchos_ParameterList.hpp"
52 #include "ROL_Types.hpp"
53 
54 namespace ROL {
55 
56 template<class Real>
57 struct SecantState {
58  std::vector<Teuchos::RCP<Vector<Real> > > iterDiff; // Step Storage
59  std::vector<Teuchos::RCP<Vector<Real> > > gradDiff; // Gradient Storage
60  std::vector<Real> product; // Step-Gradient Inner Product Storage
61  std::vector<Real> product2; // Step-Gradient Inner Product Storage
62  int storage; // Storage Size
63  int current; // Current Storage Size
64  int iter; // Current Optimization Iteration
65 };
66 
67 template<class Real>
68 class Secant {
69 private:
70 
71  Teuchos::RCP<SecantState<Real> > state_; // Secant State
72 
73 public:
74 
75  virtual ~Secant() {}
76 
77  // Constructor
78  Secant( int M = 10 ) {
79  state_ = Teuchos::rcp( new SecantState<Real> );
80  state_->storage = M;
81  state_->current = -1;
82  state_->iter = 0;
83  }
84 
85  Teuchos::RCP<SecantState<Real> >& get_state() { return state_; }
86 
87  // Update Secant Approximation
88  virtual void update( const Vector<Real> &grad, const Vector<Real> &gp, const Vector<Real> &s,
89  const Real snorm, const int iter ) {
90  state_->iter = iter;
91  Teuchos::RCP<Vector<Real> > gradDiff = grad.clone();
92  gradDiff->set(grad);
93  gradDiff->axpy(-1.0,gp);
94 
95  Real sy = s.dot(gradDiff->dual());
96  if (sy > ROL_EPSILON*snorm*snorm) {
97  if (state_->current < state_->storage-1) {
98  state_->current++; // Increment Storage
99  }
100  else {
101  state_->iterDiff.erase(state_->iterDiff.begin()); // Remove first element of s list
102  state_->gradDiff.erase(state_->gradDiff.begin()); // Remove first element of y list
103  state_->product.erase(state_->product.begin()); // Remove first element of rho list
104  }
105  state_->iterDiff.push_back(s.clone());
106  state_->iterDiff[state_->current]->set(s); // s=x_{k+1}-x_k
107  state_->gradDiff.push_back(grad.clone());
108  state_->gradDiff[state_->current]->set(*gradDiff); // y=g_{k+1}-g_k
109  state_->product.push_back(sy); // ys=1/rho
110  }
111  }
112 
113  // Apply Secant Approximate Inverse Hessian
114  virtual void applyH( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x ) = 0;
115 
116  // Apply Initial Secant Approximate Inverse Hessian
117  virtual void applyH0( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x ) {
118  Hv.set(v.dual());
119  if (state_->iter != 0 && state_->current != -1) {
120  Real yy = state_->gradDiff[state_->current]->dot(*(state_->gradDiff[state_->current]));
121  Hv.scale(state_->product[state_->current]/yy);
122  }
123  }
124 
125  // Apply Secant Approximate Hessian
126  virtual void applyB( Vector<Real> &Bv, const Vector<Real> &v, const Vector<Real> &x ) = 0;
127 
128  // Apply Initial Secant Approximate Hessian
129  virtual void applyB0( Vector<Real> &Bv, const Vector<Real> &v, const Vector<Real> &x ) {
130  Bv.set(v.dual());
131  if (state_->iter != 0 && state_->current != -1) {
132  Real yy = state_->gradDiff[state_->current]->dot(*(state_->gradDiff[state_->current]));
133  Bv.scale(yy/state_->product[state_->current]);
134  }
135  }
136 
137  // Test Secant Approximations
138  void test( const Vector<Real> &x, const Vector<Real> &s ) {
139  Teuchos::RCP<Vector<Real> > vec = x.clone();
140  Teuchos::RCP<Vector<Real> > Hvec = x.clone();
141  Teuchos::RCP<Vector<Real> > Bvec = x.clone();
142 
143  // Print BHv -> Should be v
144  vec->set(s);
145  applyH(*Hvec,*vec,x);
146  applyB(*Bvec,*Hvec,x);
147  vec->axpy(-1.0,*Bvec);
148  std::cout << " ||BHv-v|| = " << vec->norm() << "\n";
149 
150  // Print HBv -> Should be v
151  vec->set(s);
152  applyB(*Bvec,*vec,x);
153  applyH(*Hvec,*Bvec,x);
154  vec->axpy(-1.0,*Hvec);
155  std::cout << " ||HBv-v|| = " << vec->norm() << "\n";
156  }
157 
158 };
159 
160 }
161 
162 #include "ROL_SecantFactory.hpp"
163 
164 #endif
virtual void update(const Vector< Real > &grad, const Vector< Real > &gp, const Vector< Real > &s, const Real snorm, const int iter)
Definition: ROL_Secant.hpp:88
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
std::vector< Teuchos::RCP< Vector< Real > > > iterDiff
Definition: ROL_Secant.hpp:58
virtual void scale(const Real alpha)=0
Compute where .
virtual void applyB0(Vector< Real > &Bv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_Secant.hpp:129
virtual void applyH0(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_Secant.hpp:117
std::vector< Teuchos::RCP< Vector< Real > > > gradDiff
Definition: ROL_Secant.hpp:59
void test(const Vector< Real > &x, const Vector< Real > &s)
Definition: ROL_Secant.hpp:138
Contains definitions of custom data types in ROL.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
virtual Real dot(const Vector &x) const =0
Compute where .
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:68
std::vector< Real > product2
Definition: ROL_Secant.hpp:61
Secant(int M=10)
Definition: ROL_Secant.hpp:78
Teuchos::RCP< SecantState< Real > > & get_state()
Definition: ROL_Secant.hpp:85
virtual ~Secant()
Definition: ROL_Secant.hpp:75
Teuchos::RCP< SecantState< Real > > state_
Definition: ROL_Secant.hpp:71
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
std::vector< Real > product
Definition: ROL_Secant.hpp:60
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118