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
54namespace ROL {
55
56template<class Real>
57class lSR1 : public Secant<Real> {
58private:
59
60 //mutable bool updateIterate_;
62 mutable bool H0called_, B0called_;
63 Ptr<Vector<Real>> Bs_, Hy_, prim_, dual_;
64
65 using Secant<Real>::state_;
66 using Secant<Real>::y_;
68 using Secant<Real>::Bscaling_;
69
70public:
71 lSR1(int M, bool useDefaultScaling = true, Real Bscaling = Real(1), ESecantMode mode = SECANTMODE_BOTH)
72 : Secant<Real>(M,useDefaultScaling,Bscaling,mode), isInitialized_(false),
73 H0called_(false), B0called_(false) {
74 if (useDefaultScaling_) Bscaling_ = static_cast<Real>(1);
75 //updateIterate_ = true;
76 }
77
78 // Update Secant Approximation
79 void updateStorage( const Vector<Real> &x, const Vector<Real> &grad,
80 const Vector<Real> &gp, const Vector<Real> &s,
81 const Real snorm, const int iter ) {
82 const Real one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
83 if ( !isInitialized_ ) {
84 state_->iterate = x.clone();
85 y_ = grad.clone();
86 if (state_->mode == SECANTMODE_FORWARD) {
87 Bs_ = grad.clone(); dual_ = grad.clone();
88 }
89 else if (state_->mode == SECANTMODE_INVERSE) {
90 Hy_ = x.clone(); prim_ = x.clone();
91 }
92 else {
93 Bs_ = grad.clone(); dual_ = grad.clone();
94 Hy_ = x.clone(); prim_ = x.clone();
95 }
96 isInitialized_ = true;
97 }
98
99 // Update iterate
100 state_->iter = iter;
101 state_->iterate->set(x);
102
103 // Compute gradient difference
104 y_->set(grad);
105 y_->axpy(-one,gp);
106
107 Real dotF(ROL_INF<Real>()), tolF(0), dotI(ROL_INF<Real>()), tolI(0);
108 if (state_->mode == SECANTMODE_FORWARD || state_->mode == SECANTMODE_BOTH) {
109 // Compute y - Bs and <s, y - Bs>
110 applyB(*Bs_,s);
111 Bs_->scale(-one);
112 Bs_->plus(*y_);
113 //dotF = s.dot(Bs_->dual());
114 dotF = s.apply(*Bs_);
115 tolF = tol*snorm*Bs_->norm();
116 }
117 if (state_->mode == SECANTMODE_INVERSE || state_->mode == SECANTMODE_BOTH) {
118 // Compute s - Hy and <y, s - Hy>
119 applyH(*Hy_,*y_);
120 Hy_->scale(-one);
121 Hy_->plus(s);
122 //dotI = y_->dot(Hy_->dual());
123 dotI = y_->apply(*Hy_);
124 tolI = tol*y_->norm()*Hy_->norm();
125 }
126 if (std::abs(dotF) > tolF && std::abs(dotI) > tolI) {
127 if (state_->current < state_->storage-1) {
128 state_->current++;
129 if (state_->mode == SECANTMODE_INVERSE || state_->mode == SECANTMODE_BOTH) {
130 state_->iterDiff.push_back(x.clone()); // Create new memory
131 }
132 if (state_->mode == SECANTMODE_FORWARD || state_->mode == SECANTMODE_BOTH) {
133 state_->gradDiff.push_back(grad.clone()); // Create new memory
134 }
135 }
136 else {
137 if (state_->mode == SECANTMODE_INVERSE || state_->mode == SECANTMODE_BOTH) {
138 state_->iterDiff.push_back(state_->iterDiff[0]); // Move first element to the last
139 state_->iterDiff.erase(state_->iterDiff.begin()); // Remove first element of s list
140 state_->product2.erase(state_->product2.begin()); // Remove first element of rho list
141 }
142 if (state_->mode == SECANTMODE_FORWARD || state_->mode == SECANTMODE_BOTH) {
143 state_->gradDiff.push_back(state_->gradDiff[0]); // Move first element to the last
144 state_->gradDiff.erase(state_->gradDiff.begin()); // Remove first element of y list
145 state_->product.erase(state_->product.begin()); // Remove first element of rho list
146 }
147 }
148 if (state_->mode == SECANTMODE_INVERSE || state_->mode == SECANTMODE_BOTH) {
149 state_->iterDiff[state_->current]->set(*Hy_); // s_k - H_k y_k
150 state_->product2.push_back(dotI); // (s_k - H_k y_k)' y_k
151 }
152 if (state_->mode == SECANTMODE_FORWARD || state_->mode == SECANTMODE_BOTH) {
153 state_->gradDiff[state_->current]->set(*Bs_); // y_k - B_k s_k
154 state_->product.push_back(dotF); // (y_k - B_k s_k)' s_k
155 }
156 //if (useDefaultScaling_) Bscaling_ = s.dot(y_->dual())/(snorm*snorm);
157 if (useDefaultScaling_) Bscaling_ = s.apply(*y_)/(snorm*snorm);
158 }
159 /*
160 const Real one(1);
161 if ( !isInitialized_ ) {
162 state_->iterate = x.clone();
163 y_ = grad.clone();
164 isInitialized_ = true;
165 }
166
167 state_->iterate->set(x);
168 state_->iter = iter;
169 y_->set(grad);
170 y_->axpy(-one,gp);
171
172 if (updateIterate_ || state_->current == -1) {
173 Real sy = s.dot(y_->dual());
174 if (state_->current < state_->storage-1) {
175 state_->current++; // Increment Storage
176 state_->iterDiff.push_back(s.clone()); // Create new memory
177 state_->gradDiff.push_back(grad.clone()); // Create new memory
178 }
179 else {
180 state_->iterDiff.push_back(state_->iterDiff[0]); // Move first element to the last
181 state_->gradDiff.push_back(state_->gradDiff[0]); // Move first element to the last
182 state_->iterDiff.erase(state_->iterDiff.begin()); // Remove first element of s list
183 state_->gradDiff.erase(state_->gradDiff.begin()); // Remove first element of y list
184 state_->product.erase(state_->product.begin()); // Remove first element of rho list
185 }
186 state_->iterDiff[state_->current]->set(s); // s=x_{k+1}-x_k
187 state_->gradDiff[state_->current]->set(*y_); // y=g_{k+1}-g_k
188 state_->product.push_back(sy); // ys=1/rho
189 }
190 updateIterate_ = true;
191 */
192 }
193
194 // Apply Initial Secant Approximate Inverse Hessian
195 virtual void applyH0( Vector<Real> &Hv, const Vector<Real> &v ) const {
196 if (state_->current > -1) {
197 prim_->set(v.dual());
198 Hv.set(*prim_);
199 H0called_ = true;
200 }
201 else {
202 Hv.set(v.dual());
203 }
204 Hv.scale(static_cast<Real>(1)/Bscaling_);
205 }
206
207 // Apply lSR1 Approximate Inverse Hessian
208 void applyH( Vector<Real> &Hv, const Vector<Real> &v ) const {
209 if (state_->mode == SECANTMODE_INVERSE || state_->mode == SECANTMODE_BOTH) {
210 // Apply initial Hessian approximation to v
211 H0called_ = false;
212 applyH0(Hv,v);
213 // Apply rank one updates
214 if (state_->current > -1) {
215 Real prod(0);
216 if (!H0called_) prim_->set(v.dual());
217 for (int i = 0; i <= state_->current; ++i) {
218 prod = state_->iterDiff[i]->dot(*prim_);
219 Hv.axpy(prod/state_->product2[i],*state_->iterDiff[i]);
220 }
221 }
222 }
223 else {
224 throw Exception::NotImplemented(">>> ROL::lSR1::applyH : Not supported in forward mode!");
225 }
226 /*
227 std::vector<Ptr<Vector<Real>>> a(state_->current+1);
228 std::vector<Ptr<Vector<Real>>> b(state_->current+1);
229 Real byi(0), byj(0), bv(0), normbi(0), normyi(0), one(1);
230 for (int i = 0; i <= state_->current; i++) {
231 // Compute Hy
232 a[i] = Hv.clone();
233 applyH0(*(a[i]),*(state_->gradDiff[i]));
234 for (int j = 0; j < i; j++) {
235 byj = b[j]->dot((state_->gradDiff[j])->dual());
236 byi = b[j]->dot((state_->gradDiff[i])->dual());
237 a[i]->axpy(byi/byj,*(b[j]));
238 }
239 // Compute s - Hy
240 b[i] = Hv.clone();
241 b[i]->set(*(state_->iterDiff[i]));
242 b[i]->axpy(-one,*(a[i]));
243
244 // Compute Hv
245 byi = b[i]->dot((state_->gradDiff[i])->dual());
246 normbi = b[i]->norm();
247 normyi = (state_->gradDiff[i])->norm();
248 if ( i == state_->current && std::abs(byi) < sqrt(ROL_EPSILON<Real>())*normbi*normyi ) {
249 updateIterate_ = false;
250 }
251 else {
252 updateIterate_ = true;
253 bv = b[i]->dot(v.dual());
254 Hv.axpy(bv/byi,*(b[i]));
255 }
256 }
257 */
258 }
259
260 // Apply Initial Secant Approximate Hessian
261 virtual void applyB0( Vector<Real> &Bv, const Vector<Real> &v ) const {
262 if (state_->current > -1) {
263 dual_->set(v.dual());
264 Bv.set(*dual_);
265 B0called_ = true;
266 }
267 else {
268 Bv.set(v.dual());
269 }
270 Bv.scale(Bscaling_);
271 }
272
273 // Apply lSR1 Approximate Hessian
274 void applyB( Vector<Real> &Bv, const Vector<Real> &v ) const {
275 if (state_->mode == SECANTMODE_FORWARD || state_->mode == SECANTMODE_BOTH) {
276 // Apply initial Hessian approximation to v
277 B0called_ = false;
278 applyB0(Bv,v);
279 // Apply rank one updates
280 if (state_->current > -1) {
281 Real prod(0);
282 if (!B0called_) dual_->set(v.dual());
283 for (int i = 0; i <= state_->current; ++i) {
284 prod = state_->gradDiff[i]->dot(*dual_);
285 Bv.axpy(prod/state_->product[i],*state_->gradDiff[i]);
286 }
287 }
288 }
289 else {
290 throw Exception::NotImplemented(">>> ROL::lSR1::applyB : Not supported in inverse mode!");
291 }
292 /*
293 std::vector<Ptr<Vector<Real>>> a(state_->current+1);
294 std::vector<Ptr<Vector<Real>>> b(state_->current+1);
295 Real bsi(0), bsj(0), bv(0), normbi(0), normsi(0), one(1);
296 for (int i = 0; i <= state_->current; i++) {
297 // Compute Hy
298 a[i] = Bv.clone();
299 applyB0(*(a[i]),*(state_->iterDiff[i]));
300 for (int j = 0; j < i; j++) {
301 bsj = (state_->iterDiff[j])->dot(b[j]->dual());
302 bsi = (state_->iterDiff[i])->dot(b[j]->dual());
303 a[i]->axpy(bsi/bsj,*(b[j]));
304 }
305 // Compute s - Hy
306 b[i] = Bv.clone();
307 b[i]->set(*(state_->gradDiff[i]));
308 b[i]->axpy(-one,*(a[i]));
309
310 // Compute Hv
311 bsi = (state_->iterDiff[i])->dot(b[i]->dual());
312 normbi = b[i]->norm();
313 normsi = (state_->iterDiff[i])->norm();
314 if ( i == state_->current && std::abs(bsi) < sqrt(ROL_EPSILON<Real>())*normbi*normsi ) {
315 updateIterate_ = false;
316 }
317 else {
318 updateIterate_ = true;
319 bv = b[i]->dot(v.dual());
320 Bv.axpy(bv/bsi,*(b[i]));
321 }
322 }
323 */
324 }
325};
326
327}
328
329#endif
Contains definitions of custom data types in ROL.
Provides interface for and implements limited-memory secant operators.
bool useDefaultScaling_
const Ptr< SecantState< Real > > state_
Ptr< Vector< Real > > y_
Defines the linear algebra or vector space interface.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
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 .
Provides definitions for limited-memory SR1 operators.
Definition ROL_lSR1.hpp:57
void applyB(Vector< Real > &Bv, const Vector< Real > &v) const
Definition ROL_lSR1.hpp:274
Ptr< Vector< Real > > Bs_
Definition ROL_lSR1.hpp:63
Ptr< Vector< Real > > Hy_
Definition ROL_lSR1.hpp:63
Ptr< Vector< Real > > dual_
Definition ROL_lSR1.hpp:63
bool isInitialized_
Definition ROL_lSR1.hpp:61
void applyH(Vector< Real > &Hv, const Vector< Real > &v) const
Definition ROL_lSR1.hpp:208
Ptr< Vector< Real > > prim_
Definition ROL_lSR1.hpp:63
bool B0called_
Definition ROL_lSR1.hpp:62
virtual void applyH0(Vector< Real > &Hv, const Vector< Real > &v) const
Definition ROL_lSR1.hpp:195
lSR1(int M, bool useDefaultScaling=true, Real Bscaling=Real(1), ESecantMode mode=SECANTMODE_BOTH)
Definition ROL_lSR1.hpp:71
void updateStorage(const Vector< Real > &x, const Vector< Real > &grad, const Vector< Real > &gp, const Vector< Real > &s, const Real snorm, const int iter)
Definition ROL_lSR1.hpp:79
virtual void applyB0(Vector< Real > &Bv, const Vector< Real > &v) const
Definition ROL_lSR1.hpp:261
bool H0called_
Definition ROL_lSR1.hpp:62
ESecantMode
@ SECANTMODE_FORWARD
@ SECANTMODE_INVERSE
@ SECANTMODE_BOTH