ROL
ROL_TypeB_NewtonKrylovAlgorithm_Def.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_TYPEB_NEWTONKRYLOVALGORITHM_DEF_HPP
45#define ROL_TYPEB_NEWTONKRYLOVALGORITHM_DEF_HPP
46
47namespace ROL {
48namespace TypeB {
49
50template<typename Real>
52 const Ptr<Secant<Real>> &secant)
53 : secant_(secant), esec_(SECANT_USERDEFINED) {
55
56 if ( secant_ == nullPtr ) {
57 secantName_ = list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS");
59 secant_ = SecantFactory<Real>(list);
60 }
61 else {
62 secantName_ = list.sublist("General").sublist("Secant").get("User Defined Secant Name",
63 "Unspecified User Defined Secant Method");
64 }
65
66 krylovName_ = list.sublist("General").sublist("Krylov").get("Type","Conjugate Gradients");
68 krylov_ = KrylovFactory<Real>(list);
69}
70
71template<typename Real>
73 const Ptr<Krylov<Real>> &krylov,
74 const Ptr<Secant<Real>> &secant)
75 : secant_(secant), esec_(SECANT_USERDEFINED), krylov_(krylov), ekv_(KRYLOV_USERDEFINED) {
77
78 if ( secant_ == nullPtr ) {
79 secantName_ = list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS");
81 secant_ = SecantFactory<Real>(list);
82 }
83 else {
84 secantName_ = list.sublist("General").sublist("Secant").get("User Defined Secant Name",
85 "Unspecified User Defined Secant Method");
86 }
87
88 krylovName_ = list.sublist("General").sublist("Krylov").get("User Defined Krylov Name",
89 "Unspecified User Defined Krylov Method");
90}
91
92template<typename Real>
94 // Set status test
95 status_->reset();
96 status_->add(makePtr<StatusTest<Real>>(list));
97
98 // Parse parameter list
99 ParameterList &lslist = list.sublist("Step").sublist("Line Search");
100 maxit_ = lslist.get("Function Evaluation Limit", 20);
101 alpha0_ = lslist.get("Initial Step Size", 1.0);
102 useralpha_ = lslist.get("User Defined Initial Step Size", false);
103 usePrevAlpha_ = lslist.get("Use Previous Step Length as Initial Guess", false);
104 c1_ = lslist.get("Sufficient Decrease Tolerance", 1e-4);
105 rhodec_ = lslist.sublist("Line-Search Method").get("Backtracking Rate", 0.5);
106
107 useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
108 useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
109
110 verbosity_ = list.sublist("General").get("Output Level", 0);
111 writeHeader_ = verbosity_ > 2;
112}
113
114template<typename Real>
116 const Vector<Real> &g,
117 Objective<Real> &obj,
119 std::ostream &outStream) {
120 const Real one(1);
121 if (proj_ == nullPtr) {
122 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
123 }
124 // Initialize data
126 // Update approximate gradient and approximate objective function.
127 Real ftol = std::sqrt(ROL_EPSILON<Real>());
128 proj_->project(x,outStream);
129 state_->iterateVec->set(x);
130 obj.update(x,UpdateType::Initial,state_->iter);
131 state_->value = obj.value(x,ftol); state_->nfval++;
132 obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
133 state_->stepVec->set(x);
134 state_->stepVec->axpy(-one,state_->gradientVec->dual());
135 proj_->project(*state_->stepVec,outStream);
136 state_->stepVec->axpy(-one,x);
137 state_->gnorm = state_->stepVec->norm();
138 state_->snorm = ROL_INF<Real>();
139 if (!useralpha_) alpha0_ = one;
140 state_->searchSize = alpha0_;
141}
142
143template<typename Real>
145 const Vector<Real> &g,
146 Objective<Real> &obj,
148 std::ostream &outStream ) {
149 const Real one(1);
150 // Initialize trust-region data
151 initialize(x,g,obj,bnd,outStream);
152 Ptr<Vector<Real>> s = x.clone(), gp = x.clone(), gold = g.clone();
153 Ptr<Vector<Real>> pwa = x.clone(), pwa1 = x.clone();
154 Real ftrial(0), gs(0), tol(std::sqrt(ROL_EPSILON<Real>()));
155
156 Ptr<LinearOperator<Real>> hessian, precond;
157
158 // Output
159 if (verbosity_ > 0) writeOutput(outStream,true);
160
161 // Compute steepest descent step
162 gp->set(state_->gradientVec->dual());
163 while (status_->check(*state_)) {
164 // Compute step
165 hessian = makePtr<HessianPNK>(makePtrFromRef(obj),makePtrFromRef(bnd),
166 state_->iterateVec,state_->gradientVec,state_->gnorm,
167 secant_,useSecantHessVec_,pwa);
168 precond = makePtr<PrecondPNK>(makePtrFromRef(obj),makePtrFromRef(bnd),
169 state_->iterateVec,state_->gradientVec,state_->gnorm,
170 secant_,useSecantPrecond_,pwa1);
171 flagKrylov_ = 0;
172 krylov_->run(*s,*hessian,*state_->gradientVec,*precond,iterKrylov_,flagKrylov_);
173 if (flagKrylov_ == 2 && iterKrylov_ <= 1) {
174 s->set(*gp);
175 }
176 // Perform backtracking line search
177 if (!usePrevAlpha_) state_->searchSize = alpha0_;
178 x.set(*state_->iterateVec);
179 x.axpy(-state_->searchSize,*s);
180 proj_->project(x,outStream);
182 ftrial = obj.value(x,tol); ls_nfval_ = 1;
183 state_->stepVec->set(x);
184 state_->stepVec->axpy(-one,*state_->iterateVec);
185 gs = state_->stepVec->dot(*gp);
186 if (verbosity_ > 1) {
187 outStream << " In TypeB::NewtonKrylovAlgorithm: Line Search" << std::endl;
188 outStream << " Step size: " << state_->searchSize << std::endl;
189 outStream << " Trial objective value: " << ftrial << std::endl;
190 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
191 outStream << " Dot product of gradient and step: " << gs << std::endl;
192 outStream << " Sufficient decrease bound: " << -gs*c1_ << std::endl;
193 outStream << " Number of function evaluations: " << ls_nfval_ << std::endl;
194 }
195 while ( state_->value - ftrial < -c1_*gs && ls_nfval_ < maxit_ ) {
196 state_->searchSize *= rhodec_;
197 x.set(*state_->iterateVec);
198 x.axpy(-state_->searchSize,*s);
199 proj_->project(x,outStream);
201 ftrial = obj.value(x,tol); ls_nfval_++;
202 state_->stepVec->set(x);
203 state_->stepVec->axpy(-one,*state_->iterateVec);
204 gs = state_->stepVec->dot(*gp);
205 if (verbosity_ > 1) {
206 outStream << std::endl;
207 outStream << " Step size: " << state_->searchSize << std::endl;
208 outStream << " Trial objective value: " << ftrial << std::endl;
209 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
210 outStream << " Dot product of gradient and step: " << gs << std::endl;
211 outStream << " Sufficient decrease bound: " << -gs*c1_ << std::endl;
212 outStream << " Number of function evaluations: " << ls_nfval_ << std::endl;
213 }
214 }
215 state_->nfval += ls_nfval_;
216
217 // Compute norm of step
218 state_->snorm = state_->stepVec->norm();
219
220 // Update iterate
221 state_->iterateVec->set(x);
222
223 // Compute new value and gradient
224 state_->iter++;
225 state_->value = ftrial;
226 obj.update(x,UpdateType::Accept,state_->iter);
227 gold->set(*state_->gradientVec);
228 obj.gradient(*state_->gradientVec,x,tol); state_->ngrad++;
229 gp->set(state_->gradientVec->dual());
230
231 // Compute projected gradient norm
232 s->set(x); s->axpy(-one,*gp);
233 proj_->project(*s,outStream);
234 s->axpy(-one,x);
235 state_->gnorm = s->norm();
236
237 // Update secant
238 secant_->updateStorage(x,*state_->gradientVec,*gold,*state_->stepVec,state_->snorm,state_->iter);
239
240 // Update Output
241 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
242 }
243 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
244}
245
246template<typename Real>
248 std::ostream &outStream ) {
249 if (problem.getPolyhedralProjection() == nullPtr) {
250 return TypeB::Algorithm<Real>::run(problem,outStream);
251 }
252 else {
253 throw Exception::NotImplemented(">>> TypeB::NewtonKrylovAlgorithm::run : This algorithm cannot solve problems with linear equality constraints!");
254 }
255}
256
257template<typename Real>
259 const Vector<Real> &g,
260 Objective<Real> &obj,
262 Constraint<Real> &linear_econ,
263 Vector<Real> &linear_emul,
264 const Vector<Real> &linear_eres,
265 std::ostream &outStream ) {
266 throw Exception::NotImplemented(">>> TypeB::NewtonKrylovAlgorithm::run : This algorithm cannot solve problems with linear equality constraints!");
267}
268
269template<typename Real>
270void NewtonKrylovAlgorithm<Real>::writeHeader( std::ostream& os ) const {
271 std::stringstream hist;
272 if (verbosity_ > 1) {
273 hist << std::string(114,'-') << std::endl;
274 if (!useSecantHessVec_) {
275 hist << "Line-Search Projected Newton";
276 }
277 else {
278 hist << "Line-Search Projected Quasi-Newton with " << secantName_ << " Hessian approximation";
279 }
280 hist << " status output definitions" << std::endl << std::endl;
281 hist << " iter - Number of iterates (steps taken)" << std::endl;
282 hist << " value - Objective function value" << std::endl;
283 hist << " gnorm - Norm of the gradient" << std::endl;
284 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
285 hist << " alpha - Line search step length" << std::endl;
286 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
287 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
288 hist << " ls_#fval - Number of the times the objective function was evaluated during the line search" << std::endl;
289 hist << " iterCG - Number of Krylov iterations" << std::endl << std::endl;
290 hist << " flagGC - Krylov flag" << std::endl;
291 for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
292 hist << " " << NumberToString(flag) << " - "
293 << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
294 }
295 hist << std::string(114,'-') << std::endl;
296 }
297
298 hist << " ";
299 hist << std::setw(6) << std::left << "iter";
300 hist << std::setw(15) << std::left << "value";
301 hist << std::setw(15) << std::left << "gnorm";
302 hist << std::setw(15) << std::left << "snorm";
303 hist << std::setw(15) << std::left << "alpha";
304 hist << std::setw(10) << std::left << "#fval";
305 hist << std::setw(10) << std::left << "#grad";
306 hist << std::setw(10) << std::left << "#ls_fval";
307 hist << std::setw(10) << std::left << "iterCG";
308 hist << std::setw(10) << std::left << "flagCG";
309 hist << std::endl;
310 os << hist.str();
311}
312
313template<typename Real>
314void NewtonKrylovAlgorithm<Real>::writeName( std::ostream& os ) const {
315 std::stringstream hist;
316 if (!useSecantHessVec_) {
317 hist << std::endl << "Line-Search Projected Newton (Type B, Bound Constraints)" << std::endl;
318 }
319 else {
320 hist << std::endl << "Line-Search Projected Quasi-Newton with "
321 << secantName_ << " Hessian approximation" << std::endl;
322 }
323 os << hist.str();
324}
325
326template<typename Real>
327void NewtonKrylovAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
328 std::stringstream hist;
329 hist << std::scientific << std::setprecision(6);
330 if ( state_->iter == 0 ) writeName(os);
331 if ( write_header ) writeHeader(os);
332 if ( state_->iter == 0 ) {
333 hist << " ";
334 hist << std::setw(6) << std::left << state_->iter;
335 hist << std::setw(15) << std::left << state_->value;
336 hist << std::setw(15) << std::left << state_->gnorm;
337 hist << std::setw(15) << std::left << "---";
338 hist << std::setw(15) << std::left << "---";
339 hist << std::setw(10) << std::left << state_->nfval;
340 hist << std::setw(10) << std::left << state_->ngrad;
341 hist << std::setw(10) << std::left << "---";
342 hist << std::setw(10) << std::left << "---";
343 hist << std::setw(10) << std::left << "---";
344 hist << std::endl;
345 }
346 else {
347 hist << " ";
348 hist << std::setw(6) << std::left << state_->iter;
349 hist << std::setw(15) << std::left << state_->value;
350 hist << std::setw(15) << std::left << state_->gnorm;
351 hist << std::setw(15) << std::left << state_->snorm;
352 hist << std::setw(15) << std::left << state_->searchSize;
353 hist << std::setw(10) << std::left << state_->nfval;
354 hist << std::setw(10) << std::left << state_->ngrad;
355 hist << std::setw(10) << std::left << ls_nfval_;
356 hist << std::setw(10) << std::left << iterKrylov_;
357 hist << std::setw(10) << std::left << flagKrylov_;
358 hist << std::endl;
359 }
360 os << hist.str();
361}
362
363} // namespace TypeB
364} // namespace ROL
365
366#endif
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:58
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
Provides an interface to check status of optimization algorithms.
virtual void run(Problem< Real > &problem, std::ostream &outStream=std::cout)
Run algorithm on bound constrained problems (Type-B). This is the primary Type-B interface.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void writeExitStatus(std::ostream &os) const
Ptr< Secant< Real > > secant_
Secant object (used for quasi-Newton)
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout) override
Run algorithm on bound constrained problems (Type-B). This general interface supports the use of dual...
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
NewtonKrylovAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
void writeName(std::ostream &os) const override
Print step name.
Ptr< Krylov< Real > > krylov_
Krylov solver object (used for inexact Newton)
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void writeHeader(std::ostream &os) const override
Print iterate header.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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 .
Definition: ROL_Vector.hpp:153
EKrylov StringToEKrylov(std::string s)
std::string NumberToString(T Number)
Definition: ROL_Types.hpp:81
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:543
@ SECANT_USERDEFINED
Definition: ROL_Types.hpp:491
@ CG_FLAG_UNDEFINED
Definition: ROL_Types.hpp:827
@ CG_FLAG_SUCCESS
Definition: ROL_Types.hpp:822
@ KRYLOV_USERDEFINED
std::string ECGFlagToString(ECGFlag cgf)
Definition: ROL_Types.hpp:831