ROL
ROL_TypeB_QuasiNewtonAlgorithm_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_QUASINEWTONALGORITHM_DEF_HPP
45#define ROL_TYPEB_QUASINEWTONALGORITHM_DEF_HPP
46
53#include "ROL_PQNObjective.hpp"
54
55namespace ROL {
56namespace TypeB {
57
58template<typename Real>
60 const Ptr<Secant<Real>> &secant)
61 : secant_(secant), esec_(SECANT_USERDEFINED), list_(list), hasLEC_(true) {
62 // Set status test
63 status_->reset();
64 status_->add(makePtr<StatusTest<Real>>(list));
65
66 // Parse parameter list
67 ParameterList &lslist = list.sublist("Step").sublist("Line Search");
68 maxit_ = lslist.get("Function Evaluation Limit", 20);
69 c1_ = lslist.get("Sufficient Decrease Tolerance", 1e-4);
70 rhodec_ = lslist.sublist("Line-Search Method").get("Backtracking Rate", 0.5);
71 sigma1_ = lslist.sublist("PQN").get("Lower Step Size Safeguard", 0.1);
72 sigma2_ = lslist.sublist("PQN").get("Upper Step Size Safeguard", 0.9);
73 algoName_ = lslist.sublist("PQN").get("Subproblem Solver","Spectral Gradient");
74 int sp_maxit = lslist.sublist("PQN").get("Subproblem Iteration Limit", 1000);
75 sp_tol1_ = lslist.sublist("PQN").get("Subproblem Absolute Tolerance", 1e-4);
76 sp_tol2_ = lslist.sublist("PQN").get("Subproblem Relative Tolerance", 1e-2);
77 Real opt_tol = lslist.sublist("Status Test").get("Gradient Tolerance", 1e-8);
78 sp_tol_min_ = static_cast<Real>(1e-2)*opt_tol;
79 verbosity_ = list.sublist("General").get("Output Level", 0);
81
82 list_.sublist("Status Test").set("Iteration Limit", sp_maxit);
83 list_.sublist("General").set("Output Level", verbosity_>0 ? verbosity_-1 : 0);
84
85 if ( secant_ == nullPtr ) {
86 secantName_ = list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS");
88 secant_ = SecantFactory<Real>(list);
89 }
90 else {
91 secantName_ = list.sublist("General").sublist("Secant").get("User Defined Secant Name",
92 "Unspecified User Defined Secant Method");
93 }
94}
95
96
97template<typename Real>
99 const Vector<Real> &g,
100 Objective<Real> &obj,
102 std::ostream &outStream) {
103 const Real one(1);
104 if (proj_ == nullPtr) {
105 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
106 hasLEC_ = false;
107 }
108 // Initialize data
110 // Update approximate gradient and approximate objective function.
111 Real ftol = std::sqrt(ROL_EPSILON<Real>());
112 proj_->project(x,outStream); state_->nproj++;
113 state_->iterateVec->set(x);
114 obj.update(x,UpdateType::Initial,state_->iter);
115 state_->value = obj.value(x,ftol); state_->nfval++;
116 obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
117 state_->stepVec->set(x);
118 state_->stepVec->axpy(-one,state_->gradientVec->dual());
119 proj_->project(*state_->stepVec,outStream); state_->nproj++;
120 state_->stepVec->axpy(-one,x);
121 state_->gnorm = state_->stepVec->norm();
122 state_->snorm = ROL_INF<Real>();
123}
124
125template<typename Real>
127 const Vector<Real> &g,
128 Objective<Real> &obj,
130 std::ostream &outStream ) {
131 const Real half(0.5), one(1);
132 // Initialize trust-region data
133 initialize(x,g,obj,bnd,outStream);
134 Ptr<Vector<Real>> s = x.clone(), gp = x.clone(), gold = g.clone(), xs = x.clone();
135 Real ftrial(0), gs(0), alphaTmp(0), tol(std::sqrt(ROL_EPSILON<Real>())), gtol(1);
136
137 Ptr<TypeB::Algorithm<Real>> algo;
138 Ptr<PQNObjective<Real>> qobj = makePtr<PQNObjective<Real>>(secant_,x,g);
139 Ptr<Problem<Real>> problem = makePtr<Problem<Real>>(qobj,xs);
140 problem->addBoundConstraint(makePtrFromRef(bnd));
141 if (hasLEC_) {
142 problem->addLinearConstraint("LEC",proj_->getLinearConstraint(),
143 proj_->getMultiplier(),
144 proj_->getResidual());
145 problem->setProjectionAlgorithm(list_);
146 }
147 problem->finalize(false,verbosity_>2,outStream);
148
149 // Output
150 if (verbosity_ > 0) writeOutput(outStream,true);
151
152 // Compute steepest descent step
153 gp->set(state_->gradientVec->dual());
154 while (status_->check(*state_)) {
155 // Compute step
156 qobj->setAnchor(x,*state_->gradientVec);
157 xs->set(x); xs->axpy(-one,*gp); proj_->project(*xs,outStream); state_->nproj++;
158 gtol = std::max(sp_tol_min_,std::min(sp_tol1_,sp_tol2_*state_->gnorm));
159 list_.sublist("Status Test").set("Gradient Tolerance",gtol);
160 if (algoName_ == "Trust Region") algo = makePtr<TypeB::LinMoreAlgorithm<Real>>(list_);
161 else if (algoName_ == "Line Search") algo = makePtr<TypeB::GradientAlgorithm<Real>>(list_);
162 else if (algoName_ == "Primal Dual Active Set") algo = makePtr<TypeB::PrimalDualActiveSetAlgorithm<Real>>(list_);
163 else if (algoName_ == "Moreau-Yosida") algo = makePtr<TypeB::MoreauYosidaAlgorithm<Real>>(list_);
164 else if (algoName_ == "Interior Point") algo = makePtr<TypeB::InteriorPointAlgorithm<Real>>(list_);
165 else algo = makePtr<TypeB::SpectralGradientAlgorithm<Real>>(list_);
166 algo->run(*problem,outStream);
167 s->set(*xs); s->axpy(-one,x);
168 spgIter_ = algo->getState()->iter;
169 state_->nproj += staticPtrCast<const TypeB::AlgorithmState<Real>>(algo->getState())->nproj;
170
171 // Perform backtracking line search
172 state_->searchSize = one;
173 x.set(*state_->iterateVec);
174 x.axpy(state_->searchSize,*s);
176 ftrial = obj.value(x,tol); ls_nfval_ = 1;
177 gs = state_->gradientVec->apply(*s);
178 if (verbosity_ > 1) {
179 outStream << " In TypeB::QuasiNewtonAlgorithm: Line Search" << std::endl;
180 outStream << " Step size: " << state_->searchSize << std::endl;
181 outStream << " Trial objective value: " << ftrial << std::endl;
182 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
183 outStream << " Dot product of gradient and step: " << gs << std::endl;
184 outStream << " Sufficient decrease bound: " << -gs*state_->searchSize*c1_ << std::endl;
185 outStream << " Number of function evaluations: " << ls_nfval_ << std::endl;
186 }
187 while ( ftrial > state_->value + c1_*state_->searchSize*gs && ls_nfval_ < maxit_ ) {
188 alphaTmp = -half*state_->searchSize*state_->searchSize*gs
189 / (ftrial-state_->value-state_->searchSize*gs);
190 state_->searchSize = (sigma1_*state_->searchSize <= alphaTmp && alphaTmp <= sigma2_*state_->searchSize)
191 ? alphaTmp : rhodec_*state_->searchSize;
192 //state_->searchSize *= rhodec_;
193 x.set(*state_->iterateVec);
194 x.axpy(state_->searchSize,*s);
196 ftrial = obj.value(x,tol); ls_nfval_++;
197 if (verbosity_ > 1) {
198 outStream << std::endl;
199 outStream << " Step size: " << state_->searchSize << std::endl;
200 outStream << " Trial objective value: " << ftrial << std::endl;
201 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
202 outStream << " Dot product of gradient and step: " << gs << std::endl;
203 outStream << " Sufficient decrease bound: " << -gs*state_->searchSize*c1_ << std::endl;
204 outStream << " Number of function evaluations: " << ls_nfval_ << std::endl;
205 }
206 }
207 state_->nfval += ls_nfval_;
208
209 // Compute norm of step
210 state_->stepVec->set(*s);
211 state_->stepVec->scale(state_->searchSize);
212 state_->snorm = state_->stepVec->norm();
213
214 // Update iterate
215 state_->iterateVec->set(x);
216
217 // Compute new value and gradient
218 state_->iter++;
219 state_->value = ftrial;
220 obj.update(x,UpdateType::Accept,state_->iter);
221 gold->set(*state_->gradientVec);
222 obj.gradient(*state_->gradientVec,x,tol); state_->ngrad++;
223 gp->set(state_->gradientVec->dual());
224
225 // Compute projected gradient norm
226 s->set(x); s->axpy(-one,*gp);
227 proj_->project(*s,outStream); state_->nproj++;
228 s->axpy(-one,x);
229 state_->gnorm = s->norm();
230
231 // Update secant
232 secant_->updateStorage(x,*state_->gradientVec,*gold,*state_->stepVec,state_->snorm,state_->iter);
233
234 // Update Output
235 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
236 }
237 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
238}
239
240template<typename Real>
241void QuasiNewtonAlgorithm<Real>::writeHeader( std::ostream& os ) const {
242 std::stringstream hist;
243 if (verbosity_ > 1) {
244 hist << std::string(114,'-') << std::endl;
245 hist << "Line-Search Projected Quasi-Newton with " << secantName_ << " Hessian approximation";
246 hist << " status output definitions" << std::endl << std::endl;
247 hist << " iter - Number of iterates (steps taken)" << std::endl;
248 hist << " value - Objective function value" << std::endl;
249 hist << " gnorm - Norm of the gradient" << std::endl;
250 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
251 hist << " alpha - Line search step length" << std::endl;
252 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
253 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
254 hist << " #proj - Cumulative number of times the projection was computed" << std::endl;
255 hist << " ls_#fval - Number of the times the objective function was evaluated during the line search" << std::endl;
256 hist << " sp_iter - Number iterations to compute quasi-Newton step" << std::endl;
257 hist << std::string(114,'-') << std::endl;
258 }
259
260 hist << " ";
261 hist << std::setw(6) << std::left << "iter";
262 hist << std::setw(15) << std::left << "value";
263 hist << std::setw(15) << std::left << "gnorm";
264 hist << std::setw(15) << std::left << "snorm";
265 hist << std::setw(15) << std::left << "alpha";
266 hist << std::setw(10) << std::left << "#fval";
267 hist << std::setw(10) << std::left << "#grad";
268 hist << std::setw(10) << std::left << "#proj";
269 hist << std::setw(10) << std::left << "#ls_fval";
270 hist << std::setw(10) << std::left << "sp_iter";
271 hist << std::endl;
272 os << hist.str();
273}
274
275template<typename Real>
276void QuasiNewtonAlgorithm<Real>::writeName( std::ostream& os ) const {
277 std::stringstream hist;
278 hist << std::endl << "Line-Search Projected Quasi-Newton (Type B, Bound Constraints)" << std::endl;
279 os << hist.str();
280}
281
282template<typename Real>
283void QuasiNewtonAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
284 std::stringstream hist;
285 hist << std::scientific << std::setprecision(6);
286 if ( state_->iter == 0 ) writeName(os);
287 if ( write_header ) writeHeader(os);
288 if ( state_->iter == 0 ) {
289 hist << " ";
290 hist << std::setw(6) << std::left << state_->iter;
291 hist << std::setw(15) << std::left << state_->value;
292 hist << std::setw(15) << std::left << state_->gnorm;
293 hist << std::setw(15) << std::left << "---";
294 hist << std::setw(15) << std::left << "---";
295 hist << std::setw(10) << std::left << state_->nfval;
296 hist << std::setw(10) << std::left << state_->ngrad;
297 hist << std::setw(10) << std::left << state_->nproj;
298 hist << std::setw(10) << std::left << "---";
299 hist << std::setw(10) << std::left << "---";
300 hist << std::endl;
301 }
302 else {
303 hist << " ";
304 hist << std::setw(6) << std::left << state_->iter;
305 hist << std::setw(15) << std::left << state_->value;
306 hist << std::setw(15) << std::left << state_->gnorm;
307 hist << std::setw(15) << std::left << state_->snorm;
308 hist << std::setw(15) << std::left << state_->searchSize;
309 hist << std::setw(10) << std::left << state_->nfval;
310 hist << std::setw(10) << std::left << state_->ngrad;
311 hist << std::setw(10) << std::left << state_->nproj;
312 hist << std::setw(10) << std::left << ls_nfval_;
313 hist << std::setw(10) << std::left << spgIter_;
314 hist << std::endl;
315 }
316 os << hist.str();
317}
318
319} // namespace TypeB
320} // namespace ROL
321
322#endif
Provides the interface to apply upper and lower bound constraints.
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.
Provides interface for and implements limited-memory secant operators.
Provides an interface to check status of optimization algorithms.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void writeExitStatus(std::ostream &os) const
const Ptr< CombinedStatusTest< Real > > status_
Real sigma2_
Upper safeguard for quadratic line search (default: 0.9)
Real c1_
Sufficient Decrease Parameter (default: 1e-4)
int maxit_
Maximum number of line search steps (default: 20)
Real sigma1_
Lower safeguard for quadratic line search (default: 0.1)
Ptr< Secant< Real > > secant_
Secant object (used for quasi-Newton)
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
void writeHeader(std::ostream &os) const override
Print iterate header.
void writeName(std::ostream &os) const override
Print step name.
Real rhodec_
Backtracking rate (default: 0.5)
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...
QuasiNewtonAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
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 .
ESecant StringToESecant(std::string s)
@ SECANT_USERDEFINED