ROL
ROL_TypeU_BundleAlgorithm_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 parlist of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this parlist 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_TYPEU_BUNDLEALGORITHM_DEF_H
45 #define ROL_TYPEU_BUNDLEALGORITHM_DEF_H
46 
47 #include "ROL_BundleStatusTest.hpp"
48 #include "ROL_Bundle_U_AS.hpp"
49 #include "ROL_Bundle_U_TT.hpp"
51 
52 namespace ROL {
53 namespace TypeU {
54 
55 
56 template<typename Real>
58  const Ptr<LineSearch_U<Real>> &lineSearch )
59  : Algorithm<Real>(),
60  bundle_(ROL::nullPtr), lineSearch_(ROL::nullPtr),
61  QPiter_(0), QPmaxit_(0), QPtol_(0), step_flag_(0),
62  T_(0), tol_(0), m1_(0), m2_(0), m3_(0), nu_(0),
63  ls_maxit_(0), first_print_(true), isConvex_(false) {
64  // Set bundle status test
65  status_->reset();
66  status_->add(makePtr<BundleStatusTest<Real>>(parlist));
67 
68  // Parse parameter parlist
69  const Real zero(0), two(2), oem3(1.e-3), oem6(1.e-6), oem8(1.e-8);
70  const Real p1(0.1), p2(0.2), p9(0.9), oe3(1.e3), oe8(1.e8);
71  ParameterList &blist = parlist.sublist("Step").sublist("Bundle");
72  state_->searchSize = blist.get("Initial Trust-Region Parameter", oe3);
73  T_ = blist.get("Maximum Trust-Region Parameter", oe8);
74  tol_ = blist.get("Epsilon Solution Tolerance", oem6);
75  m1_ = blist.get("Upper Threshold for Serious Step", p1);
76  m2_ = blist.get("Lower Threshold for Serious Step", p2);
77  m3_ = blist.get("Upper Threshold for Null Step", p9);
78  nu_ = blist.get("Tolerance for Trust-Region Parameter", oem3);
79 
80  // Initialize bundle
81  Real coeff = blist.get("Distance Measure Coefficient", zero);
82  Real omega = blist.get("Locality Measure Coefficient", two);
83  unsigned maxSize = blist.get("Maximum Bundle Size", 200);
84  unsigned remSize = blist.get("Removal Size for Bundle Update", 2);
85  if ( blist.get("Cutting Plane Solver",0) == 1 ) {
86  bundle_ = makePtr<Bundle_U_TT<Real>>(maxSize,coeff,omega,remSize);
87  }
88  else {
89  bundle_ = makePtr<Bundle_U_AS<Real>>(maxSize,coeff,omega,remSize);
90  }
91  isConvex_ = ((coeff == zero) ? true : false);
92 
93  // Initialize QP solver
94  QPtol_ = blist.get("Cutting Plane Tolerance", oem8);
95  QPmaxit_ = blist.get("Cutting Plane Iteration Limit", 1000);
96 
97  // Initialize Line Search
98  ParameterList &lslist = parlist.sublist("Step").sublist("Line Search");
99  ls_maxit_ = lslist.get("Maximum Number of Function Evaluations",20);
100  if ( !isConvex_ && lineSearch_==nullPtr ) {
101  lineSearch_ = LineSearchUFactory<Real>(parlist);
102  }
103 
104  // Get verbosity level
105  verbosity_ = parlist.sublist("General").get("Output Level", 0);
106  printHeader_ = verbosity_ > 2;
107 }
108 
109 template<typename Real>
111  const Vector<Real> &g,
112  Objective<Real> &obj,
113  std::ostream &outStream) {
114  // Initialize data
116  if (!isConvex_) {
117  lineSearch_->initialize(x,g);
118  }
119  // Update objective function, get value and gradient
120  Real tol = std::sqrt(ROL_EPSILON<Real>());
121  obj.update(x,UpdateType::Initial,state_->iter);
122  state_->value = obj.value(x,tol);
123  state_->nfval++;
124  obj.gradient(*state_->gradientVec,x,tol);
125  state_->ngrad++;
126  state_->gnorm = state_->gradientVec->norm();
127  bundle_->initialize(*state_->gradientVec);
128 }
129 
130 template<typename Real>
132  const Vector<Real> &g,
133  Objective<Real> &obj,
134  std::ostream &outStream ) {
135  const Real zero(0), two(2), half(0.5);
136  // Initialize trust-region data
137  Real tol(std::sqrt(ROL_EPSILON<Real>()));
138  initialize(x,g,obj,outStream);
139  Ptr<Vector<Real>> y = x.clone();
140  Ptr<Vector<Real>> aggSubGradNew = g.clone();
141  Real aggSubGradOldNorm = state_->gnorm;
142  Real linErrNew(0), valueNew(0);
143  Real aggLinErrNew(0), aggLinErrOld(0), aggDistMeasNew(0);
144  Real v(0), l(0), u(T_), gd(0);
145  bool flag = true;
146 
147  // Output
148  if (verbosity_ > 0) writeOutput(outStream,true);
149 
150  while (status_->check(*state_)) {
151  first_print_ = false; // Print header only on first serious step
152  QPiter_ = (step_flag_==1 ? 0 : QPiter_); // Reset QPiter only on serious steps
153  v = zero; l = zero; u = T_; gd = zero;
154  flag = true;
155  while (flag) {
156  /*************************************************************/
157  /******** Solve Dual Cutting Plane QP Problem ****************/
158  /*************************************************************/
159  QPiter_ += bundle_->solveDual(state_->searchSize,QPmaxit_,QPtol_); // Solve QP subproblem
160  bundle_->aggregate(*aggSubGradNew,aggLinErrNew,aggDistMeasNew); // Compute aggregate info
161  state_->aggregateGradientNorm = aggSubGradNew->norm(); // Aggregate subgradient norm
162  if (verbosity_ > 1) {
163  outStream << std::endl;
164  outStream << " Computation of aggregrate quantities" << std::endl;
165  outStream << " Aggregate subgradient norm: " << state_->aggregateGradientNorm << std::endl;
166  outStream << " Aggregate linearization error: " << aggLinErrNew << std::endl;
167  outStream << " Aggregate distance measure: " << aggDistMeasNew << std::endl;
168  }
169  /*************************************************************/
170  /******** Construct Cutting Plane Solution *******************/
171  /*************************************************************/
172  v = -state_->searchSize*std::pow(state_->aggregateGradientNorm,two)-aggLinErrNew; // CP objective value
173  state_->stepVec->set(aggSubGradNew->dual());
174  state_->stepVec->scale(-state_->searchSize); // CP solution
175  state_->snorm = state_->searchSize*state_->aggregateGradientNorm; // Step norm
176  if (verbosity_ > 1) {
177  outStream << std::endl;
178  outStream << " Solve cutting plan subproblem" << std::endl;
179  outStream << " Cutting plan objective value: " << v << std::endl;
180  outStream << " Norm of computed step: " << state_->snorm << std::endl;
181  outStream << " 'Trust-region' radius: " << state_->searchSize << std::endl;
182  }
183  /*************************************************************/
184  /******** Decide Whether Step is Serious or Null *************/
185  /*************************************************************/
186  if (std::max(state_->aggregateGradientNorm,aggLinErrNew) <= tol_) {
187  // Current iterate is already epsilon optimal!
188  state_->stepVec->zero(); state_->snorm = zero;
189  flag = false;
190  step_flag_ = 1;
191  state_->flag = true;
192  break;
193  }
194  else if (std::isnan(state_->aggregateGradientNorm)
195  || std::isnan(aggLinErrNew)
196  || (std::isnan(aggDistMeasNew) && !isConvex_)) {
197  state_->stepVec->zero(); state_->snorm = zero;
198  flag = false;
199  step_flag_ = 2;
200  state_->flag = true;
201  }
202  else {
203  // Current iterate is not epsilon optimal.
204  y->set(x); y->plus(*state_->stepVec); // y is the candidate iterate
205  obj.update(*y,UpdateType::Accept,state_->iter); // Update objective at y
206  valueNew = obj.value(*y,tol); // Compute objective value at y
207  state_->nfval++;
208  obj.gradient(*state_->gradientVec,*y,tol); // Compute objective (sub)gradient at y
209  state_->ngrad++;
210  // Compute new linearization error and distance measure
211  //gd = state_->stepVec->dot(state_->gradientVec->dual());
212  gd = state_->stepVec->apply(*state_->gradientVec);
213  linErrNew = state_->value - (valueNew - gd); // Linearization error
214  // Determine whether to take a serious or null step
215  Real eps = static_cast<Real>(10)*ROL_EPSILON<Real>();
216  Real del = eps*std::max(static_cast<Real>(1),std::abs(state_->value));
217  Real Df = (valueNew - state_->value) - del;
218  Real Dm = v - del;
219  bool SS1 = false;
220  if (std::abs(Df) < eps && std::abs(Dm) < eps) {
221  SS1 = true;
222  }
223  else {
224  SS1 = (Df < m1_*Dm);
225  }
226  //bool SS1 = (valueNew-state_->value < m1_*v);
227  //bool NS1 = (valueNew-state_->value >= m1_*v);
228  bool NS2a = (bundle_->computeAlpha(state_->snorm,linErrNew) <= m3_*aggLinErrOld);
229  bool NS2b = (std::abs(state_->value-valueNew) <= aggSubGradOldNorm + aggLinErrOld);
230  if (verbosity_ > 1) {
231  outStream << std::endl;
232  outStream << " Check for serious/null step" << std::endl;
233  outStream << " Serious step test SS(i): " << SS1 << std::endl;
234  outStream << " -> Left hand side: " << valueNew-state_->value << std::endl;
235  outStream << " -> Right hand side: " << m1_*v << std::endl;
236  outStream << " Null step test NS(iia): " << NS2a << std::endl;
237  outStream << " -> Left hand side: " << bundle_->computeAlpha(state_->snorm,linErrNew) << std::endl;
238  outStream << " -> Right hand side: " << m3_*aggLinErrOld << std::endl;
239  outStream << " Null step test NS(iib): " << NS2b << std::endl;
240  outStream << " -> Left hand side: " << std::abs(state_->value-valueNew) << std::endl;
241  outStream << " -> Right hand side: " << aggSubGradOldNorm + aggLinErrOld << std::endl;
242  }
243  if ( isConvex_ ) {
244  /************* Convex objective ****************/
245  if ( SS1 ) {
246  bool SS2 = (gd >= m2_*v || state_->searchSize >= T_-nu_);
247  if (verbosity_ > 1) {
248  outStream << " Serious step test SS(iia): " << (gd >= m2_*v) << std::endl;
249  outStream << " -> Left hand side: " << gd << std::endl;
250  outStream << " -> Right hand side: " << m2_*v << std::endl;
251  outStream << " Serious step test SS(iia): " << (state_->searchSize >= T_-nu_) << std::endl;
252  outStream << " -> Left hand side: " << state_->searchSize << std::endl;
253  outStream << " -> Right hand side: " << T_-nu_ << std::endl;
254  }
255  if ( SS2 ) { // Serious Step
256  step_flag_ = 1;
257  flag = false;
258  if (verbosity_ > 1) {
259  outStream << " Serious step taken" << std::endl;
260  }
261  break;
262  }
263  else { // Increase trust-region radius
264  l = state_->searchSize;
265  state_->searchSize = half*(u+l);
266  if (verbosity_ > 1) {
267  outStream << " Increase 'trust-region' radius: " << state_->searchSize << std::endl;
268  }
269  }
270  }
271  else {
272  if ( NS2a || NS2b ) { // Null step
273  state_->stepVec->zero(); state_->snorm = zero;
274  step_flag_ = 0;
275  flag = false;
276  if (verbosity_ > 1) {
277  outStream << " Null step taken" << std::endl;
278  }
279  break;
280  }
281  else { // Decrease trust-region radius
282  u = state_->searchSize;
283  state_->searchSize = half*(u+l);
284  if (verbosity_ > 1) {
285  outStream << " Decrease 'trust-region' radius: " << state_->searchSize << std::endl;
286  }
287  }
288  }
289  }
290  else {
291  /************* Nonconvex objective *************/
292  bool NS3 = (gd - bundle_->computeAlpha(state_->snorm,linErrNew) >= m2_*v);
293  if (verbosity_ > 1) {
294  outStream << " Null step test NS(iii): " << NS3 << std::endl;
295  outStream << " -> Left hand side: " << gd - bundle_->computeAlpha(state_->snorm,linErrNew) << std::endl;
296  outStream << " -> Right hand side: " << m2_*v << std::endl;
297  }
298  if ( SS1 ) { // Serious step
299  step_flag_ = 1;
300  flag = false;
301  break;
302  }
303  else {
304  if ( NS2a || NS2b ) {
305  if ( NS3 ) { // Null step
306  state_->stepVec->zero();
307  step_flag_ = 0;
308  flag = false;
309  break;
310  }
311  else {
312  if ( NS2b ) { // Line search
313  Real alpha = zero;
314  int ls_nfval = 0, ls_ngrad = 0;
315  lineSearch_->run(alpha,valueNew,ls_nfval,ls_ngrad,gd,*state_->stepVec,x,obj);
316  if ( ls_nfval == ls_maxit_ ) { // Null step
317  state_->stepVec->zero();
318  step_flag_ = 0;
319  flag = false;
320  break;
321  }
322  else { // Serious step
323  state_->stepVec->scale(alpha);
324  step_flag_ = 1;
325  flag = false;
326  break;
327  }
328  }
329  else { // Decrease trust-region radius
330  u = state_->searchSize;
331  state_->searchSize = half*(u+l);
332  }
333  }
334  }
335  else { // Decrease trust-region radius
336  u = state_->searchSize;
337  state_->searchSize = half*(u+l);
338  }
339  }
340  }
341  }
342  } // End While
343  /*************************************************************/
344  /******** Update Algorithm State *****************************/
345  /*************************************************************/
346  state_->aggregateModelError = aggLinErrNew;
347  aggSubGradOldNorm = state_->aggregateGradientNorm;
348  aggLinErrOld = aggLinErrNew;
349 
350  if ( !state_->flag ) {
351  /*************************************************************/
352  /******** Reset Bundle If Maximum Size Reached ***************/
353  /*************************************************************/
354  bundle_->reset(*aggSubGradNew,aggLinErrNew,state_->snorm);
355  /*************************************************************/
356  /******** Update Bundle with Step Information ****************/
357  /*************************************************************/
358  if ( step_flag_==1 ) {
359  // Serious step was taken
360  x.plus(*state_->stepVec); // Update iterate
361  Real valueOld = state_->value; // Store previous objective value
362  state_->value = valueNew; // Add new objective value to state
363  bundle_->update(step_flag_,valueNew-valueOld,state_->snorm,*state_->gradientVec,*state_->stepVec);
364  }
365  else if ( step_flag_==0 ) {
366  // Null step was taken
367  bundle_->update(step_flag_,linErrNew,state_->snorm,*state_->gradientVec,*state_->stepVec);
368  }
369  }
370  /*************************************************************/
371  /******** Update Algorithm State *****************************/
372  /*************************************************************/
373  state_->iterateVec->set(x);
374  state_->gnorm = state_->gradientVec->norm();
375  if ( step_flag_==1 ) {
376  state_->iter++;
377  }
378 
379  // Update Output
380  if (verbosity_ > 0) writeOutput(outStream,printHeader_);
381  }
382  if (verbosity_ > 0) Algorithm<Real>::writeExitStatus(outStream);
383 }
384 
385 template<typename Real>
386 void BundleAlgorithm<Real>::writeHeader( std::ostream& os ) const {
387  std::ios_base::fmtflags osFlags(os.flags());
388  os << " ";
389  os << std::setw(6) << std::left << "iter";
390  os << std::setw(15) << std::left << "value";
391  os << std::setw(15) << std::left << "gnorm";
392  os << std::setw(15) << std::left << "snorm";
393  os << std::setw(10) << std::left << "#fval";
394  os << std::setw(10) << std::left << "#grad";
395  os << std::setw(15) << std::left << "znorm";
396  os << std::setw(15) << std::left << "alpha";
397  os << std::setw(15) << std::left << "TRparam";
398  os << std::setw(10) << std::left << "QPiter";
399  os << std::endl;
400  os.flags(osFlags);
401 }
402 
403 template<typename Real>
404 void BundleAlgorithm<Real>::writeName(std::ostream& os) const {
405  std::ios_base::fmtflags osFlags(os.flags());
406  os << std::endl << "Bundle Trust-Region Algorithm" << std::endl;
407  os.flags(osFlags);
408 }
409 
410 template<typename Real>
411 void BundleAlgorithm<Real>::writeOutput( std::ostream& os, bool print_header) const {
412  std::ios_base::fmtflags osFlags(os.flags());
413  os << std::scientific << std::setprecision(6);
414  if ( state_->iter == 0 && first_print_ ) {
415  writeName(os);
416  if ( print_header ) {
417  writeHeader(os);
418  }
419  os << " ";
420  os << std::setw(6) << std::left << state_->iter;
421  os << std::setw(15) << std::left << state_->value;
422  os << std::setw(15) << std::left << state_->gnorm;
423  os << std::setw(15) << std::left << "---";
424  os << std::setw(10) << std::left << state_->nfval;
425  os << std::setw(10) << std::left << state_->ngrad;
426  os << std::setw(15) << std::left << "---";
427  os << std::setw(15) << std::left << "---";
428  os << std::setw(15) << std::left << state_->searchSize;
429  os << std::setw(10) << std::left << "---";
430  os << std::endl;
431  }
432  if ( step_flag_==1 && state_->iter > 0 ) {
433  if ( print_header ) {
434  writeHeader(os);
435  }
436  else {
437  os << " ";
438  os << std::setw(6) << std::left << state_->iter;
439  os << std::setw(15) << std::left << state_->value;
440  os << std::setw(15) << std::left << state_->gnorm;
441  os << std::setw(15) << std::left << state_->snorm;
442  os << std::setw(10) << std::left << state_->nfval;
443  os << std::setw(10) << std::left << state_->ngrad;
444  os << std::setw(15) << std::left << state_->aggregateGradientNorm;
445  os << std::setw(15) << std::left << state_->aggregateModelError;
446  os << std::setw(15) << std::left << state_->searchSize;
447  os << std::setw(10) << std::left << QPiter_;
448  os << std::endl;
449  }
450  }
451  os.flags(osFlags);
452 }
453 
454 } // namespace TypeU
455 } // namespace ROL
456 
457 #endif
Provides interface for and implements line searches.
Provides the interface to evaluate objective functions.
void writeName(std::ostream &os) const override
Print step name.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
BundleAlgorithm(ParameterList &parlist, const Ptr< LineSearch_U< Real >> &lineSearch=nullPtr)
virtual void plus(const Vector &x)=0
Compute , where .
const Ptr< AlgorithmState< Real > > state_
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void initialize(const Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, std::ostream &outStream=std::cout)
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void writeHeader(std::ostream &os) const override
Print iterate header.
Provides an interface to run unconstrained optimization algorithms.
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, std::ostream &outStream=std::cout) override
Run algorithm on unconstrained problems (Type-U). This general interface supports the use of dual opt...
void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
virtual void writeExitStatus(std::ostream &os) const
const Ptr< CombinedStatusTest< Real > > status_
Ptr< LineSearch_U< Real > > lineSearch_