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