ROL
ROL_TypeB_GradientAlgorithm_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_GRADIENTALGORITHM_DEF_HPP
45#define ROL_TYPEB_GRADIENTALGORITHM_DEF_HPP
46
47namespace ROL {
48namespace TypeB {
49
50template<typename Real>
52 // Set status test
53 status_->reset();
54 status_->add(makePtr<StatusTest<Real>>(list));
55
56 // Parse parameter list
57 ParameterList &lslist = list.sublist("Step").sublist("Line Search");
58 maxit_ = lslist.get("Function Evaluation Limit", 20);
59 alpha0_ = lslist.get("Initial Step Size", 1.0);
60 normAlpha_ = lslist.get("Normalize Initial Step Size", false);
61 alpha0bnd_ = lslist.get("Lower Bound for Initial Step Size", 1e-4);
62 useralpha_ = lslist.get("User Defined Initial Step Size", false);
63 usePrevAlpha_ = lslist.get("Use Previous Step Length as Initial Guess", false);
64 c1_ = lslist.get("Sufficient Decrease Tolerance", 1e-4);
65 maxAlpha_ = lslist.get("Maximum Step Size", alpha0_);
66 useAdapt_ = lslist.get("Use Adaptive Step Size Selection", true);
67 rhodec_ = lslist.sublist("Line-Search Method").get("Backtracking Rate", 0.5);
68 rhoinc_ = lslist.sublist("Line-Search Method").get("Increase Rate" , 2.0);
69 verbosity_ = list.sublist("General").get("Output Level", 0);
70 writeHeader_ = verbosity_ > 2;
71}
72
73template<typename Real>
75 const Vector<Real> &g,
76 Objective<Real> &obj,
78 std::ostream &outStream) {
79 const Real one(1);
80 if (proj_ == nullPtr) {
81 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
82 }
83 // Initialize data
85 // Update approximate gradient and approximate objective function.
86 Real ftol = std::sqrt(ROL_EPSILON<Real>());
87 proj_->project(x,outStream);
88 obj.update(x,UpdateType::Initial,state_->iter);
89 state_->value = obj.value(x,ftol);
90 state_->nfval++;
91 obj.gradient(*state_->gradientVec,x,ftol);
92 state_->ngrad++;
93 state_->stepVec->set(x);
94 state_->stepVec->axpy(-one,state_->gradientVec->dual());
95 proj_->project(*state_->stepVec,outStream);
96 Real fnew = state_->value;
97 if (!useralpha_) {
98 // Evaluate objective at P(x - g)
99 obj.update(*state_->stepVec,UpdateType::Trial);
100 fnew = obj.value(*state_->stepVec,ftol);
102 state_->nfval++;
103 }
104 state_->stepVec->axpy(-one,x);
105 state_->gnorm = state_->stepVec->norm();
106 state_->snorm = ROL_INF<Real>();
107 if (!useralpha_) {
108 const Real half(0.5);
109 // Minimize quadratic interpolate to compute new alpha
110 //Real gs = state_->stepVec->dot(state_->gradientVec->dual());
111 Real gs = state_->stepVec->apply(*state_->gradientVec);
112 Real denom = (fnew - state_->value - gs);
113 bool flag = maxAlpha_ == alpha0_;
114 alpha0_ = ((denom > ROL_EPSILON<Real>()) ? -half*gs/denom : alpha0bnd_);
115 alpha0_ = ((alpha0_ > alpha0bnd_) ? alpha0_ : one);
116 if (flag) maxAlpha_ = alpha0_;
117 }
118 // Normalize initial CP step length
119 if (normAlpha_) {
120 alpha0_ /= state_->gradientVec->norm();
121 }
122 state_->searchSize = alpha0_;
123}
124
125template<typename Real>
127 const Vector<Real> &g,
128 Objective<Real> &obj,
130 std::ostream &outStream ) {
131 const Real one(1);
132 // Initialize trust-region data
133 initialize(x,g,obj,bnd,outStream);
134 Ptr<Vector<Real>> s = x.clone();
135 Real ftrial(0), gs(0), ftrialP(0), alphaP(0), tol(std::sqrt(ROL_EPSILON<Real>()));
136 int ls_nfval = 0;
137 bool incAlpha = false, accept = true;
138
139 // Output
140 if (verbosity_ > 0) writeOutput(outStream,true);
141
142 // Compute steepest descent step
143 state_->stepVec->set(state_->gradientVec->dual());
144 while (status_->check(*state_)) {
145 accept = true;
146 // Perform backtracking line search
147 if (!usePrevAlpha_ && !useAdapt_) state_->searchSize = alpha0_;
148 state_->iterateVec->set(x);
149 state_->iterateVec->axpy(-state_->searchSize,*state_->stepVec);
150 proj_->project(*state_->iterateVec,outStream);
151 obj.update(*state_->iterateVec,UpdateType::Trial);
152 ftrial = obj.value(*state_->iterateVec,tol);
153 ls_nfval = 1;
154 s->set(*state_->iterateVec);
155 s->axpy(-one,x);
156 gs = s->dot(*state_->stepVec);
157 incAlpha = (state_->value - ftrial >= -c1_*gs);
158 if (verbosity_ > 1) {
159 outStream << " In TypeB::GradientAlgorithm: Line Search" << std::endl;
160 outStream << " Step size: " << state_->searchSize << std::endl;
161 outStream << " Trial objective value: " << ftrial << std::endl;
162 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
163 outStream << " Dot product of gradient and step: " << gs << std::endl;
164 outStream << " Sufficient decrease bound: " << -gs*c1_ << std::endl;
165 outStream << " Number of function evaluations: " << ls_nfval << std::endl;
166 outStream << " Increase alpha?: " << incAlpha << std::endl;
167 }
168 if (incAlpha && useAdapt_) {
169 ftrialP = ROL_INF<Real>();
170 while ( state_->value - ftrial >= -c1_*gs
171 && ftrial <= ftrialP
172 && state_->searchSize < maxAlpha_
173 && ls_nfval < maxit_ ) {
174 // Previous value was acceptable
175 obj.update(*state_->iterateVec,UpdateType::Accept);
176 alphaP = state_->searchSize;
177 ftrialP = ftrial;
178 state_->searchSize *= rhoinc_;
179 state_->searchSize = std::min(state_->searchSize,maxAlpha_);
180 state_->iterateVec->set(x);
181 state_->iterateVec->axpy(-state_->searchSize,*state_->stepVec);
182 proj_->project(*state_->iterateVec,outStream);
183 obj.update(*state_->iterateVec,UpdateType::Trial);
184 ftrial = obj.value(*state_->iterateVec,tol);
185 ls_nfval++;
186 s->set(*state_->iterateVec);
187 s->axpy(-one,x);
188 gs = s->dot(*state_->stepVec);
189 if (verbosity_ > 1) {
190 outStream << std::endl;
191 outStream << " Step size: " << state_->searchSize << std::endl;
192 outStream << " Trial objective value: " << ftrial << std::endl;
193 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
194 outStream << " Dot product of gradient and step: " << gs << std::endl;
195 outStream << " Sufficient decrease bound: " << -gs*c1_ << std::endl;
196 outStream << " Number of function evaluations: " << ls_nfval << std::endl;
197 }
198 }
199 if (state_->value - ftrial < -c1_*gs || ftrial > ftrialP) {
200 ftrial = ftrialP;
201 state_->searchSize = alphaP;
202 state_->iterateVec->set(x);
203 state_->iterateVec->axpy(-state_->searchSize,*state_->stepVec);
204 proj_->project(*state_->iterateVec,outStream);
205 s->set(*state_->iterateVec);
206 s->axpy(-one,x);
207 accept = false;
208 }
209 }
210 else {
211 while ( state_->value - ftrial < -c1_*gs && ls_nfval < maxit_ ) {
212 state_->searchSize *= rhodec_;
213 state_->iterateVec->set(x);
214 state_->iterateVec->axpy(-state_->searchSize,*state_->stepVec);
215 proj_->project(*state_->iterateVec,outStream);
216 obj.update(*state_->iterateVec,UpdateType::Trial);
217 ftrial = obj.value(*state_->iterateVec,tol);
218 ls_nfval++;
219 s->set(*state_->iterateVec);
220 s->axpy(-one,x);
221 gs = s->dot(*state_->stepVec);
222 if (verbosity_ > 1) {
223 outStream << std::endl;
224 outStream << " Step size: " << state_->searchSize << std::endl;
225 outStream << " Trial objective value: " << ftrial << std::endl;
226 outStream << " Computed reduction: " << state_->value-ftrial << std::endl;
227 outStream << " Dot product of gradient and step: " << gs << std::endl;
228 outStream << " Sufficient decrease bound: " << -gs*c1_ << std::endl;
229 outStream << " Number of function evaluations: " << ls_nfval << std::endl;
230 }
231 }
232 }
233 state_->nfval += ls_nfval;
234
235 // Compute norm of step
236 state_->stepVec->set(*s);
237 state_->snorm = state_->stepVec->norm();
238
239 // Update iterate
240 x.set(*state_->iterateVec);
241
242 // Compute new value and gradient
243 state_->iter++;
244 state_->value = ftrial;
245 if (accept) obj.update(x,UpdateType::Accept,state_->iter);
246 else obj.update(x,UpdateType::Revert,state_->iter);
247 obj.gradient(*state_->gradientVec,x,tol);
248 state_->ngrad++;
249
250 // Compute steepest descent step
251 state_->stepVec->set(state_->gradientVec->dual());
252
253 // Compute projected gradient norm
254 s->set(x); s->axpy(-one,*state_->stepVec);
255 proj_->project(*s,outStream);
256 s->axpy(-one,x);
257 state_->gnorm = s->norm();
258
259 // Update Output
260 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
261 }
262 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
263}
264
265template<typename Real>
266void GradientAlgorithm<Real>::writeHeader( std::ostream& os ) const {
267 std::stringstream hist;
268 if (verbosity_ > 1) {
269 hist << std::string(109,'-') << std::endl;
270 hist << "Projected gradient descent";
271 hist << " status output definitions" << std::endl << std::endl;
272 hist << " iter - Number of iterates (steps taken)" << std::endl;
273 hist << " value - Objective function value" << std::endl;
274 hist << " gnorm - Norm of the gradient" << std::endl;
275 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
276 hist << " alpha - Line search step length" << std::endl;
277 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
278 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
279 hist << std::string(109,'-') << std::endl;
280 }
281
282 hist << " ";
283 hist << std::setw(6) << std::left << "iter";
284 hist << std::setw(15) << std::left << "value";
285 hist << std::setw(15) << std::left << "gnorm";
286 hist << std::setw(15) << std::left << "snorm";
287 hist << std::setw(15) << std::left << "alpha";
288 hist << std::setw(10) << std::left << "#fval";
289 hist << std::setw(10) << std::left << "#grad";
290 hist << std::endl;
291 os << hist.str();
292}
293
294template<typename Real>
295void GradientAlgorithm<Real>::writeName( std::ostream& os ) const {
296 std::stringstream hist;
297 hist << std::endl << "Projected Gradient Descent with Backtracking Line Search (Type B, Bound Constraints)" << std::endl;
298 os << hist.str();
299}
300
301template<typename Real>
302void GradientAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
303 std::stringstream hist;
304 hist << std::scientific << std::setprecision(6);
305 if ( state_->iter == 0 ) writeName(os);
306 if ( write_header ) writeHeader(os);
307 if ( state_->iter == 0 ) {
308 hist << " ";
309 hist << std::setw(6) << std::left << state_->iter;
310 hist << std::setw(15) << std::left << state_->value;
311 hist << std::setw(15) << std::left << state_->gnorm;
312 hist << std::setw(15) << std::left << "---";
313 hist << std::setw(15) << std::left << "---";
314 hist << std::setw(10) << std::left << state_->nfval;
315 hist << std::setw(10) << std::left << state_->ngrad;
316 hist << std::endl;
317 }
318 else {
319 hist << " ";
320 hist << std::setw(6) << std::left << state_->iter;
321 hist << std::setw(15) << std::left << state_->value;
322 hist << std::setw(15) << std::left << state_->gnorm;
323 hist << std::setw(15) << std::left << state_->snorm;
324 hist << std::setw(15) << std::left << state_->searchSize;
325 hist << std::setw(10) << std::left << state_->nfval;
326 hist << std::setw(10) << std::left << state_->ngrad;
327 hist << std::endl;
328 }
329 os << hist.str();
330}
331
332} // namespace TypeB
333} // namespace ROL
334
335#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 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
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 writeHeader(std::ostream &os) const override
Print iterate header.
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void writeName(std::ostream &os) const override
Print step name.
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
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.