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.
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.
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.
virtual void set(const Vector &x)
Set where .
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 .
EKrylov StringToEKrylov(std::string s)
std::string NumberToString(T Number)
Definition ROL_Types.hpp:81
ESecant StringToESecant(std::string s)
@ SECANT_USERDEFINED
@ CG_FLAG_UNDEFINED
@ CG_FLAG_SUCCESS
std::string ECGFlagToString(ECGFlag cgf)