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 
47 namespace ROL {
48 namespace TypeB {
49 
50 template<typename Real>
51 GradientAlgorithm<Real>::GradientAlgorithm(ParameterList &list) {
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 
73 template<typename Real>
74 void GradientAlgorithm<Real>::initialize(Vector<Real> &x,
75  const Vector<Real> &g,
76  Objective<Real> &obj,
77  BoundConstraint<Real> &bnd,
78  std::ostream &outStream) {
79  const Real one(1);
80  if (proj_ == nullPtr) {
81  proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
82  }
83  // Initialize data
84  TypeB::Algorithm<Real>::initialize(x,g);
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);
101  obj.update(x,UpdateType::Revert);
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 
125 template<typename Real>
126 void GradientAlgorithm<Real>::run( Vector<Real> &x,
127  const Vector<Real> &g,
128  Objective<Real> &obj,
129  BoundConstraint<Real> &bnd,
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 
265 template<typename Real>
266 void GradientAlgorithm<Real>::writeHeader( std::ostream& os ) const {
267  std::ios_base::fmtflags osFlags(os.flags());
268  if (verbosity_ > 1) {
269  os << std::string(109,'-') << std::endl;
270  os << "Projected gradient descent";
271  os << " status output definitions" << std::endl << std::endl;
272  os << " iter - Number of iterates (steps taken)" << std::endl;
273  os << " value - Objective function value" << std::endl;
274  os << " gnorm - Norm of the gradient" << std::endl;
275  os << " snorm - Norm of the step (update to optimization vector)" << std::endl;
276  os << " alpha - Line search step length" << std::endl;
277  os << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
278  os << " #grad - Cumulative number of times the gradient was computed" << std::endl;
279  os << std::string(109,'-') << std::endl;
280  }
281 
282  os << " ";
283  os << std::setw(6) << std::left << "iter";
284  os << std::setw(15) << std::left << "value";
285  os << std::setw(15) << std::left << "gnorm";
286  os << std::setw(15) << std::left << "snorm";
287  os << std::setw(15) << std::left << "alpha";
288  os << std::setw(10) << std::left << "#fval";
289  os << std::setw(10) << std::left << "#grad";
290  os << std::endl;
291  os.flags(osFlags);
292 }
293 
294 template<typename Real>
295 void GradientAlgorithm<Real>::writeName( std::ostream& os ) const {
296  std::ios_base::fmtflags osFlags(os.flags());
297  os << std::endl << "Projected Gradient Descent with Backtracking Line Search (Type B, Bound Constraints)" << std::endl;
298  os.flags(osFlags);
299 }
300 
301 template<typename Real>
302 void GradientAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
303  std::ios_base::fmtflags osFlags(os.flags());
304  os << std::scientific << std::setprecision(6);
305  if ( state_->iter == 0 ) writeName(os);
306  if ( write_header ) writeHeader(os);
307  if ( state_->iter == 0 ) {
308  os << " ";
309  os << std::setw(6) << std::left << state_->iter;
310  os << std::setw(15) << std::left << state_->value;
311  os << std::setw(15) << std::left << state_->gnorm;
312  os << std::setw(15) << std::left << "---";
313  os << std::setw(15) << std::left << "---";
314  os << std::setw(10) << std::left << state_->nfval;
315  os << std::setw(10) << std::left << state_->ngrad;
316  os << std::endl;
317  }
318  else {
319  os << " ";
320  os << std::setw(6) << std::left << state_->iter;
321  os << std::setw(15) << std::left << state_->value;
322  os << std::setw(15) << std::left << state_->gnorm;
323  os << std::setw(15) << std::left << state_->snorm;
324  os << std::setw(15) << std::left << state_->searchSize;
325  os << std::setw(10) << std::left << state_->nfval;
326  os << std::setw(10) << std::left << state_->ngrad;
327  os << std::endl;
328  }
329  os.flags(osFlags);
330 }
331 
332 } // namespace TypeB
333 } // namespace ROL
334 
335 #endif
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:78
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::TypeB::GradientAlgorithm::writeName
void writeName(std::ostream &os) const override
Print step name.
Definition: ROL_TypeB_GradientAlgorithm_Def.hpp:295
ROL::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::TypeB::Algorithm::writeExitStatus
virtual void writeExitStatus(std::ostream &os) const
Definition: ROL_TypeB_Algorithm_Def.hpp:304
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
ROL::Objective::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_Objective.hpp:96
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:74
ROL::TypeB::GradientAlgorithm::writeOutput
void writeOutput(std::ostream &os, const bool write_header=false) const override
Print iterate status.
Definition: ROL_TypeB_GradientAlgorithm_Def.hpp:302
ROL::TypeB::GradientAlgorithm::initialize
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
Definition: ROL_TypeB_GradientAlgorithm_Def.hpp:74
ROL::StatusTest
Provides an interface to check status of optimization algorithms.
Definition: ROL_StatusTest.hpp:58
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:73
ROL::TypeB::GradientAlgorithm::writeHeader
void writeHeader(std::ostream &os) const override
Print iterate header.
Definition: ROL_TypeB_GradientAlgorithm_Def.hpp:266
ROL::TypeB::GradientAlgorithm::run
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...
Definition: ROL_TypeB_GradientAlgorithm_Def.hpp:126
ROL::TypeB::Algorithm::initialize
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Definition: ROL_TypeB_Algorithm_Def.hpp:63
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
Generated by   doxygen 1.8.5