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
48#include "ROL_Bundle_U_AS.hpp"
49#include "ROL_Bundle_U_TT.hpp"
51
52namespace ROL {
53namespace TypeU {
54
55
56template<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);
107}
108
109template<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
130template<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
385template<typename Real>
386void BundleAlgorithm<Real>::writeHeader( std::ostream& os ) const {
387 std::stringstream hist;
388 hist << " ";
389 hist << std::setw(6) << std::left << "iter";
390 hist << std::setw(15) << std::left << "value";
391 hist << std::setw(15) << std::left << "gnorm";
392 hist << std::setw(15) << std::left << "snorm";
393 hist << std::setw(10) << std::left << "#fval";
394 hist << std::setw(10) << std::left << "#grad";
395 hist << std::setw(15) << std::left << "znorm";
396 hist << std::setw(15) << std::left << "alpha";
397 hist << std::setw(15) << std::left << "TRparam";
398 hist << std::setw(10) << std::left << "QPiter";
399 hist << std::endl;
400 os << hist.str();
401}
402
403template<typename Real>
404void BundleAlgorithm<Real>::writeName(std::ostream& os) const {
405 std::stringstream hist;
406 hist << std::endl << "Bundle Trust-Region Algorithm" << std::endl;
407 os << hist.str();
408}
409
410template<typename Real>
411void BundleAlgorithm<Real>::writeOutput( std::ostream& os, bool print_header) const {
412 std::stringstream hist;
413 hist << std::scientific << std::setprecision(6);
414 if ( state_->iter == 0 && first_print_ ) {
415 writeName(os);
416 if ( print_header ) {
417 writeHeader(os);
418 }
419 hist << " ";
420 hist << std::setw(6) << std::left << state_->iter;
421 hist << std::setw(15) << std::left << state_->value;
422 hist << std::setw(15) << std::left << state_->gnorm;
423 hist << std::setw(15) << std::left << "---";
424 hist << std::setw(10) << std::left << state_->nfval;
425 hist << std::setw(10) << std::left << state_->ngrad;
426 hist << std::setw(15) << std::left << "---";
427 hist << std::setw(15) << std::left << "---";
428 hist << std::setw(15) << std::left << state_->searchSize;
429 hist << std::setw(10) << std::left << "---";
430 hist << std::endl;
431 }
432 if ( step_flag_==1 && state_->iter > 0 ) {
433 if ( print_header ) {
434 writeHeader(os);
435 }
436 else {
437 hist << " ";
438 hist << std::setw(6) << std::left << state_->iter;
439 hist << std::setw(15) << std::left << state_->value;
440 hist << std::setw(15) << std::left << state_->gnorm;
441 hist << std::setw(15) << std::left << state_->snorm;
442 hist << std::setw(10) << std::left << state_->nfval;
443 hist << std::setw(10) << std::left << state_->ngrad;
444 hist << std::setw(15) << std::left << state_->aggregateGradientNorm;
445 hist << std::setw(15) << std::left << state_->aggregateModelError;
446 hist << std::setw(15) << std::left << state_->searchSize;
447 hist << std::setw(10) << std::left << QPiter_;
448 hist << std::endl;
449 }
450 }
451 os << hist.str();
452}
453
454} // namespace TypeU
455} // namespace ROL
456
457#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides interface for and implements line searches.
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 run unconstrained optimization algorithms.
const Ptr< CombinedStatusTest< Real > > status_
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void writeExitStatus(std::ostream &os) const
const Ptr< AlgorithmState< Real > > state_
void writeOutput(std::ostream &os, bool print_header=false) const override
Print iterate status.
BundleAlgorithm(ParameterList &parlist, const Ptr< LineSearch_U< Real > > &lineSearch=nullPtr)
Ptr< LineSearch_U< Real > > lineSearch_
void writeHeader(std::ostream &os) const override
Print iterate header.
void initialize(const Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, std::ostream &outStream=std::cout)
void writeName(std::ostream &os) const override
Print step name.
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...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.