ROL
ROL_TypeB_SpectralGradientAlgorithm_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_SPECTRALGRADIENTALGORITHM_DEF_HPP
45#define ROL_TYPEB_SPECTRALGRADIENTALGORITHM_DEF_HPP
46
47#include <deque>
48
49namespace ROL {
50namespace TypeB {
51
52template<typename Real>
54 // Set status test
55 status_->reset();
56 status_->add(makePtr<StatusTest<Real>>(list));
57
58 // Parse parameter list
59 ParameterList &lslist = list.sublist("Step").sublist("Spectral Gradient");
60 maxit_ = lslist.get("Function Evaluation Limit", 20);
61 lambda_ = lslist.get("Initial Spectral Step Size", -1.0);
62 lambdaMin_ = lslist.get("Minimum Spectral Step Size", 1e-8);
63 lambdaMax_ = lslist.get("Maximum Spectral Step Size", 1e8);
64 sigma1_ = lslist.get("Lower Step Size Safeguard", 0.1);
65 sigma2_ = lslist.get("Upper Step Size Safeguard", 0.9);
66 rhodec_ = lslist.get("Backtracking Rate", 0.5);
67 gamma_ = lslist.get("Sufficient Decrease Tolerance", 1e-4);
68 maxSize_ = lslist.get("Maximum Storage Size", 10);
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 zero(0), one(1);
80 if (proj_ == nullPtr)
81 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
82 // Initialize data
84 // Update approximate gradient and approximate objective function.
85 Real ftol = std::sqrt(ROL_EPSILON<Real>());
86 proj_->project(x,outStream); state_->nproj++;
87 obj.update(x,UpdateType::Initial,state_->iter);
88 state_->value = obj.value(x,ftol); state_->nfval++;
89 obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
90 state_->stepVec->set(x);
91 state_->stepVec->axpy(-one,state_->gradientVec->dual());
92 proj_->project(*state_->stepVec,outStream); state_->nproj++;
93 state_->stepVec->axpy(-one,x);
94 state_->gnorm = state_->stepVec->norm();
95 state_->snorm = ROL_INF<Real>();
96 if (lambda_ <= zero && state_->gnorm != zero)
97 lambda_ = std::max(lambdaMin_,std::min(one/state_->gnorm,lambdaMax_));
98}
99
100template<typename Real>
102 const Vector<Real> &g,
103 Objective<Real> &obj,
105 std::ostream &outStream ) {
106 const Real half(0.5), one(1), eps(std::sqrt(ROL_EPSILON<Real>()));
107 // Initialize trust-region data
108 initialize(x,g,obj,bnd,outStream);
109 Ptr<Vector<Real>> s = x.clone(), y = g.clone(), xmin = x.clone();
110 Real ftrial(0), fmax(0), gs(0), alpha(1), alphaTmp(1), fmin(0);
111 Real ys(0), ss(0), tol(std::sqrt(ROL_EPSILON<Real>()));
112 int ls_nfval = 0;
113 std::deque<Real> fqueue; fqueue.push_back(state_->value);
114
115 fmin = state_->value;
116 xmin->set(x);
117
118 // Output
119 if (verbosity_ > 0) writeOutput(outStream, true);
120
121 // Iterate spectral projected gradient
122 state_->stepVec->set(state_->gradientVec->dual());
123 while (status_->check(*state_)) {
124 // Compute projected spectral step
125 state_->iterateVec->set(x);
126 state_->iterateVec->axpy(-lambda_,*state_->stepVec);
127 proj_->project(*state_->iterateVec,outStream); state_->nproj++;
128 s->set(*state_->iterateVec);
129 s->axpy(-one,x);
130
131 // Nonmonotone Linesearch
132 ls_nfval = 0;
133 obj.update(*state_->iterateVec,UpdateType::Trial);
134 ftrial = obj.value(*state_->iterateVec,tol); ls_nfval++;
135 alpha = one;
136 fmax = *std::max_element(fqueue.begin(),fqueue.end());
137 gs = state_->gradientVec->apply(*s);
138 if (verbosity_ > 1) {
139 outStream << " In TypeB::SpectralGradientAlgorithm Line Search" << std::endl;
140 outStream << " Step size: " << alpha << std::endl;
141 outStream << " Trial objective value: " << ftrial << std::endl;
142 outStream << " Max stored objective value: " << fmax << std::endl;
143 outStream << " Computed reduction: " << fmax-ftrial << std::endl;
144 outStream << " Dot product of gradient and step: " << gs << std::endl;
145 outStream << " Sufficient decrease bound: " << -gs*gamma_*alpha << std::endl;
146 outStream << " Number of function evaluations: " << ls_nfval << std::endl;
147 }
148 while (ftrial > fmax + gamma_*alpha*gs && ls_nfval < maxit_) {
149 alphaTmp = -half*alpha*alpha*gs/(ftrial-state_->value-alpha*gs);
150 alpha = (sigma1_*alpha <= alphaTmp && alphaTmp <= sigma2_*alpha) ? alphaTmp : rhodec_*alpha;
151 state_->iterateVec->set(x);
152 state_->iterateVec->axpy(alpha,*s);
153 obj.update(*state_->iterateVec,UpdateType::Trial);
154 ftrial = obj.value(*state_->iterateVec,tol); ls_nfval++;
155 if (verbosity_ > 1) {
156 outStream << " In TypeB::SpectralGradientAlgorithm: Line Search" << std::endl;
157 outStream << " Step size: " << alpha << std::endl;
158 outStream << " Trial objective value: " << ftrial << std::endl;
159 outStream << " Max stored objective value: " << fmax << std::endl;
160 outStream << " Computed reduction: " << fmax-ftrial << std::endl;
161 outStream << " Dot product of gradient and step: " << gs << std::endl;
162 outStream << " Sufficient decrease bound: " << -gs*gamma_*alpha << std::endl;
163 outStream << " Number of function evaluations: " << ls_nfval << std::endl;
164 }
165 }
166 state_->nfval += ls_nfval;
167 if (static_cast<int>(fqueue.size()) == maxSize_) fqueue.pop_front();
168 fqueue.push_back(ftrial);
169
170 // Update state
171 state_->iter++;
172 state_->value = ftrial;
173 state_->searchSize = alpha;
174 x.set(*state_->iterateVec);
175 obj.update(x,UpdateType::Accept,state_->iter);
176
177 // Store the best iterate
178 if (state_->value <= fmin) {
179 fmin = state_->value;
180 xmin->set(x);
181 }
182
183 // Compute spectral step length
184 s->scale(alpha);
185 y->set(*state_->gradientVec);
186 y->scale(-one);
187 obj.gradient(*state_->gradientVec,x,tol); state_->ngrad++;
188 y->plus(*state_->gradientVec);
189 ys = y->apply(*s);
190 ss = s->dot(*s);
191 lambda_ = (ys<=eps ? lambdaMax_ : std::max(lambdaMin_,std::min(ss/ys,lambdaMax_)));
192 state_->snorm = std::sqrt(ss);
193
194 // Compute gradient step
195 state_->stepVec->set(state_->gradientVec->dual());
196
197 // Compute projected gradient norm
198 s->set(x); s->axpy(-one,*state_->stepVec);
199 proj_->project(*s,outStream); state_->nproj++;
200 s->axpy(-one,x);
201 state_->gnorm = s->norm();
202
203 // Update Output
204 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
205 }
206 x.set(*xmin);
207 state_->value = fmin;
208 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
209}
210
211template<typename Real>
212void SpectralGradientAlgorithm<Real>::writeHeader( std::ostream& os ) const {
213 std::stringstream hist;
214 if (verbosity_ > 1) {
215 hist << std::string(109,'-') << std::endl;
216 hist << "Spectral projected gradient descent";
217 hist << " status output definitions" << std::endl << std::endl;
218 hist << " iter - Number of iterates (steps taken)" << std::endl;
219 hist << " value - Objective function value" << std::endl;
220 hist << " gnorm - Norm of the gradient" << std::endl;
221 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
222 hist << " alpha - Line search step length" << std::endl;
223 hist << " lambda - Spectral step length" << std::endl;
224 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
225 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
226 hist << " #proj - Cumulative number of times the projection was computed" << std::endl;
227 hist << std::string(109,'-') << std::endl;
228 }
229
230 hist << " ";
231 hist << std::setw(6) << std::left << "iter";
232 hist << std::setw(15) << std::left << "value";
233 hist << std::setw(15) << std::left << "gnorm";
234 hist << std::setw(15) << std::left << "snorm";
235 hist << std::setw(15) << std::left << "alpha";
236 hist << std::setw(15) << std::left << "lambda";
237 hist << std::setw(10) << std::left << "#fval";
238 hist << std::setw(10) << std::left << "#grad";
239 hist << std::setw(10) << std::left << "#proj";
240 hist << std::endl;
241 os << hist.str();
242}
243
244template<typename Real>
245void SpectralGradientAlgorithm<Real>::writeName( std::ostream& os ) const {
246 std::stringstream hist;
247 hist << std::endl << "Projected Spectral Gradient Method (Type B, Bound Constraints)" << std::endl;
248 os << hist.str();
249}
250
251template<typename Real>
252void SpectralGradientAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
253 std::stringstream hist;
254 hist << std::scientific << std::setprecision(6);
255 if ( state_->iter == 0 ) writeName(os);
256 if ( write_header ) writeHeader(os);
257 if ( state_->iter == 0 ) {
258 hist << " ";
259 hist << std::setw(6) << std::left << state_->iter;
260 hist << std::setw(15) << std::left << state_->value;
261 hist << std::setw(15) << std::left << state_->gnorm;
262 hist << std::setw(15) << std::left << "---";
263 hist << std::setw(15) << std::left << "---";
264 hist << std::setw(15) << std::left << lambda_;
265 hist << std::setw(10) << std::left << state_->nfval;
266 hist << std::setw(10) << std::left << state_->ngrad;
267 hist << std::setw(10) << std::left << state_->nproj;
268 hist << std::endl;
269 }
270 else {
271 hist << " ";
272 hist << std::setw(6) << std::left << state_->iter;
273 hist << std::setw(15) << std::left << state_->value;
274 hist << std::setw(15) << std::left << state_->gnorm;
275 hist << std::setw(15) << std::left << state_->snorm;
276 hist << std::setw(15) << std::left << state_->searchSize;
277 hist << std::setw(15) << std::left << lambda_;
278 hist << std::setw(10) << std::left << state_->nfval;
279 hist << std::setw(10) << std::left << state_->ngrad;
280 hist << std::setw(10) << std::left << state_->nproj;
281 hist << std::endl;
282 }
283 os << hist.str();
284}
285
286} // namespace TypeB
287} // namespace ROL
288
289#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
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 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 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 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.