ROL
ROL_lSR1.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_LSR1_H
45 #define ROL_LSR1_H
46 
51 #include "ROL_Secant.hpp"
52 #include "ROL_Types.hpp"
53 
54 namespace ROL {
55 
56 template<class Real>
57 class lSR1 : public Secant<Real> {
58 private:
59 
61 
62 public:
63  lSR1(int M) : Secant<Real>(M) {
64  updateIterate_ = true;
65  }
66 
67  // Update Secant Approximation
68  void update( const Vector<Real> &grad, const Vector<Real> &gp, const Vector<Real> &s,
69  const Real snorm, const int iter ) {
70  // Get Generic Secant State
71  Teuchos::RCP<SecantState<Real> >& state = Secant<Real>::get_state();
72 
73  state->iter = iter;
74  Teuchos::RCP<Vector<Real> > gradDiff = grad.clone();
75  gradDiff->set(grad);
76  gradDiff->axpy(-1.0,gp);
77 
78  Real sy = s.dot(gradDiff->dual());
79  if (updateIterate_ || state->current == -1) {
80  if (state->current < state->storage-1) {
81  state->current++; // Increment Storage
82  }
83  else {
84  state->iterDiff.erase(state->iterDiff.begin()); // Remove first element of s list
85  state->gradDiff.erase(state->gradDiff.begin()); // Remove first element of y list
86  state->product.erase(state->product.begin()); // Remove first element of rho list
87  }
88  state->iterDiff.push_back(s.clone());
89  state->iterDiff[state->current]->set(s); // s=x_{k+1}-x_k
90  state->gradDiff.push_back(grad.clone());
91  state->gradDiff[state->current]->set(*gradDiff); // y=g_{k+1}-g_k
92  state->product.push_back(sy); // ys=1/rho
93  }
94  updateIterate_ = true;
95  }
96 
97  // Apply Initial Secant Approximate Inverse Hessian
98  virtual void applyH0( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x ) {
99  Hv.set(v.dual());
100  }
101 
102 
103  // Apply lSR1 Approximate Inverse Hessian
104  void applyH( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x ) {
105  // Get Generic Secant State
106  Teuchos::RCP<SecantState<Real> >& state = Secant<Real>::get_state();
107 
108  // Apply initial Hessian approximation to v
109  applyH0(Hv,v,x);
110 
111  std::vector<Teuchos::RCP<Vector<Real> > > a(state->current+1);
112  std::vector<Teuchos::RCP<Vector<Real> > > b(state->current+1);
113  Real byi = 0.0, byj = 0.0, bv = 0.0, normbi = 0.0, normyi = 0.0;
114  for (int i = 0; i <= state->current; i++) {
115  // Compute Hy
116  a[i] = Hv.clone();
117  applyH0(*(a[i]),*(state->gradDiff[i]),x);
118  for (int j = 0; j < i; j++) {
119  byj = b[j]->dot((state->gradDiff[j])->dual());
120  byi = b[j]->dot((state->gradDiff[i])->dual());
121  a[i]->axpy(byi/byj,*(b[j]));
122  }
123  // Compute s - Hy
124  b[i] = Hv.clone();
125  b[i]->set(*(state->iterDiff[i]));
126  b[i]->axpy(-1.0,*(a[i]));
127 
128  // Compute Hv
129  byi = b[i]->dot((state->gradDiff[i])->dual());
130  normbi = b[i]->norm();
131  normyi = (state->gradDiff[i])->norm();
132  if ( i == state->current && std::abs(byi) < sqrt(ROL_EPSILON)*normbi*normyi ) {
133  updateIterate_ = false;
134  }
135  else {
136  updateIterate_ = true;
137  bv = b[i]->dot(v.dual());
138  Hv.axpy(bv/byi,*(b[i]));
139  }
140  }
141  }
142 
143  // Apply Initial Secant Approximate Hessian
144  virtual void applyB0( Vector<Real> &Bv, const Vector<Real> &v, const Vector<Real> &x ) {
145  Bv.set(v.dual());
146  }
147 
148 
149  // Apply lSR1 Approximate Hessian
150  void applyB( Vector<Real> &Bv, const Vector<Real> &v, const Vector<Real> &x ) {
151  // Get Generic Secant State
152  Teuchos::RCP<SecantState<Real> >& state = Secant<Real>::get_state();
153 
154  // Apply initial Hessian approximation to v
155  applyB0(Bv,v,x);
156 
157  std::vector<Teuchos::RCP<Vector<Real> > > a(state->current+1);
158  std::vector<Teuchos::RCP<Vector<Real> > > b(state->current+1);
159  Real bsi = 0.0, bsj = 0.0, bv = 0.0, normbi = 0.0, normsi = 0.0;
160  for (int i = 0; i <= state->current; i++) {
161  // Compute Hy
162  a[i] = Bv.clone();
163  applyB0(*(a[i]),*(state->iterDiff[i]),x);
164  for (int j = 0; j < i; j++) {
165  bsj = (state->iterDiff[j])->dot(b[j]->dual());
166  bsi = (state->iterDiff[i])->dot(b[j]->dual());
167  a[i]->axpy(bsi/bsj,*(b[j]));
168  }
169  // Compute s - Hy
170  b[i] = Bv.clone();
171  b[i]->set(*(state->gradDiff[i]));
172  b[i]->axpy(-1.0,*(a[i]));
173 
174  // Compute Hv
175  bsi = (state->iterDiff[i])->dot(b[i]->dual());
176  normbi = b[i]->norm();
177  normsi = (state->iterDiff[i])->norm();
178  if ( i == state->current && std::abs(bsi) < sqrt(ROL_EPSILON)*normbi*normsi ) {
179  updateIterate_ = false;
180  }
181  else {
182  updateIterate_ = true;
183  bv = b[i]->dot(v.dual());
184  Bv.axpy(bv/bsi,*(b[i]));
185  }
186  }
187  }
188 
189 };
190 
191 }
192 
193 #endif
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
virtual void applyH0(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_lSR1.hpp:98
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
Contains definitions of custom data types in ROL.
Provides definitions for limited-memory SR1 operators.
Definition: ROL_lSR1.hpp:57
void update(const Vector< Real > &grad, const Vector< Real > &gp, const Vector< Real > &s, const Real snorm, const int iter)
Definition: ROL_lSR1.hpp:68
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
void applyB(Vector< Real > &Bv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_lSR1.hpp:150
virtual Real dot(const Vector &x) const =0
Compute where .
lSR1(int M)
Definition: ROL_lSR1.hpp:63
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:68
bool updateIterate_
Definition: ROL_lSR1.hpp:60
Teuchos::RCP< SecantState< Real > > & get_state()
Definition: ROL_Secant.hpp:85
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
virtual void applyB0(Vector< Real > &Bv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_lSR1.hpp:144
void applyH(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x)
Definition: ROL_lSR1.hpp:104
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118