ROL
ROL_TypeP_ProxGradientAlgorithm_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_TYPEP_PROXGRADIENTALGORITHM_DEF_HPP
45 #define ROL_TYPEP_PROXGRADIENTALGORITHM_DEF_HPP
46 
47 namespace ROL {
48 namespace TypeP {
49 
50 template<typename Real>
51 ProxGradientAlgorithm<Real>::ProxGradientAlgorithm(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  initProx_ = lslist.get("Apply Prox to Initial Guess", false);
68  rhodec_ = lslist.sublist("Line-Search Method").get("Backtracking Rate", 0.5);
69  rhoinc_ = lslist.sublist("Line-Search Method").get("Increase Rate" , 2.0);
70  t0_ = list.sublist("Status Test").get("Gradient Scale" , 1.0);
71  verbosity_ = list.sublist("General").get("Output Level", 0);
72  writeHeader_ = verbosity_ > 2;
73 }
74 
75 template<typename Real>
76 void ProxGradientAlgorithm<Real>::initialize(Vector<Real> &x,
77  const Vector<Real> &g,
78  Objective<Real> &sobj,
79  Objective<Real> &nobj,
80  Vector<Real> &px,
81  Vector<Real> &dg,
82  std::ostream &outStream) {
83  const Real one(1);
84  // Initialize data
85  TypeP::Algorithm<Real>::initialize(x,g);
86  // Update approximate gradient and approximate objective function.
87  Real ftol = std::sqrt(ROL_EPSILON<Real>());
88  if (initProx_) {
89  nobj.prox(*state_->iterateVec,x,state_->searchSize,ftol);
90  state_->nprox++;
91  x.set(*state_->iterateVec);
92  }
93  // Evaluate objective function
94  sobj.update(x,UpdateType::Initial,state_->iter);
95  nobj.update(x,UpdateType::Initial,state_->iter);
96  state_->svalue = sobj.value(x,ftol); state_->nsval++;
97  state_->nvalue = nobj.value(x,ftol); state_->nnval++;
98  state_->value = state_->svalue + state_->nvalue;
99  // Evaluate gradient of smooth part
100  sobj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
101  dg.set(state_->gradientVec->dual());
102  // Compute initial step size as 2/L, where L = 2|f(x+s)-f(x)-f'(x)s|/||s||^2
103  // is a lower estimate of the Lipschitz constant of f
104  if (!useralpha_) {
105  bool flag = maxAlpha_ == alpha0_;
106  // Evaluate objective at Prox(x - t0 dg)
107  pgstep(px, *state_->stepVec, nobj, x, dg, t0_, ftol);
108  state_->snorm = state_->stepVec->norm();
109  sobj.update(px,UpdateType::Trial);
110  Real snew = sobj.value(px,ftol);
111  sobj.update(x,UpdateType::Revert);
112  state_->nsval++;
113  Real gs = state_->gradientVec->apply(*state_->stepVec);
114  alpha0_ = (state_->snorm * state_->snorm) / std::abs(snew - state_->svalue - gs);
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  state_->searchSize = alpha0_;
122  // Evaluate proximal gradient
123  pgstep(*state_->iterateVec, *state_->stepVec, nobj, x, dg, state_->searchSize, ftol);
124  state_->snorm = state_->stepVec->norm();
125  state_->gnorm = state_->snorm / state_->searchSize;
126 }
127 
128 template<typename Real>
129 void ProxGradientAlgorithm<Real>::run( Vector<Real> &x,
130  const Vector<Real> &g,
131  Objective<Real> &sobj,
132  Objective<Real> &nobj,
133  std::ostream &outStream ) {
134  const Real one(1);
135  Real tol(std::sqrt(ROL_EPSILON<Real>()));
136  // Initialize trust-region data
137  Ptr<Vector<Real>> px = x.clone(), pxP = x.clone(), dg = x.clone();
138  initialize(x,g,sobj,nobj,*px,*dg,outStream);
139  Real strial(0), ntrial(0), Ftrial(0), Qk(0);
140  Real strialP(0), ntrialP(0), FtrialP(0), alphaP(0);
141  Real snorm(state_->snorm), searchSize(state_->searchSize);
142  int ls_nfval = 0;
143  bool incAlpha = false, accept = true;
144 
145  // Output
146  if (verbosity_ > 0) writeOutput(outStream,true);
147 
148  // Compute steepest descent step
149  while (status_->check(*state_)) {
150  accept = true;
151  // Perform backtracking line search
152  state_->searchSize = searchSize;
153  // Compute objective function values
154  sobj.update(*state_->iterateVec,UpdateType::Trial);
155  strial = sobj.value(*state_->iterateVec,tol);
156  nobj.update(*state_->iterateVec,UpdateType::Trial);
157  ntrial = nobj.value(*state_->iterateVec,tol);
158  Ftrial = strial + ntrial;
159  ls_nfval = 1;
160  // Compute decrease indicator
161  Qk = state_->gradientVec->apply(*state_->stepVec) + ntrial - state_->nvalue;
162  incAlpha = (Ftrial - state_->value <= c1_*Qk);
163  if (verbosity_ > 1) {
164  outStream << " In TypeP::GradientAlgorithm: Line Search" << std::endl;
165  outStream << " Step size: " << state_->searchSize << std::endl;
166  outStream << " Trial smooth value: " << strial << std::endl;
167  outStream << " Trial nonsmooth value: " << ntrial << std::endl;
168  outStream << " Computed reduction: " << state_->value-Ftrial << std::endl;
169  outStream << " Dot product of gradient and step: " << Qk << std::endl;
170  outStream << " Sufficient decrease bound: " << -Qk*c1_ << std::endl;
171  outStream << " Number of function evaluations: " << ls_nfval << std::endl;
172  outStream << " Increase alpha?: " << incAlpha << std::endl;
173  }
174  if (incAlpha && useAdapt_) {
175  ntrialP = ROL_INF<Real>();
176  strialP = ROL_INF<Real>();
177  FtrialP = ntrialP + strialP;
178  while ( Ftrial - state_->value <= c1_*Qk
179  && Ftrial <= FtrialP
180  && state_->searchSize < maxAlpha_
181  && ls_nfval < maxit_ ) {
182  // Previous value was acceptable
183  sobj.update(*state_->iterateVec,UpdateType::Accept);
184  nobj.update(*state_->iterateVec,UpdateType::Accept);
185  // Backup previous values to avoid recomputation
186  pxP->set(*state_->iterateVec);
187  alphaP = state_->searchSize;
188  strialP = strial;
189  ntrialP = ntrial;
190  FtrialP = Ftrial;
191  // Increase search size
192  state_->searchSize *= rhoinc_;
193  state_->searchSize = std::min(state_->searchSize,maxAlpha_);
194  // Compute proximal gradient step with new search size
195  pgstep(*state_->iterateVec, *state_->stepVec, nobj, x, *dg, state_->searchSize, tol);
196  // Compute objective function values
197  sobj.update(*state_->iterateVec,UpdateType::Trial);
198  strial = sobj.value(*state_->iterateVec,tol);
199  nobj.update(*state_->iterateVec,UpdateType::Trial);
200  ntrial = nobj.value(*state_->iterateVec,tol);
201  Ftrial = strial + ntrial;
202  ls_nfval++;
203  // Compute decrease indicator
204  Qk = state_->gradientVec->apply(*state_->stepVec) + ntrial - state_->nvalue;
205  if (verbosity_ > 1) {
206  outStream << std::endl;
207  outStream << " Step size: " << state_->searchSize << std::endl;
208  outStream << " Trial smooth value: " << strial << std::endl;
209  outStream << " Trial nonsmooth value: " << ntrial << std::endl;
210  outStream << " Computed reduction: " << state_->value-Ftrial << std::endl;
211  outStream << " Dot product of gradient and step: " << Qk << std::endl;
212  outStream << " Sufficient decrease bound: " << -Qk*c1_ << std::endl;
213  outStream << " Number of function evaluations: " << ls_nfval << std::endl;
214  }
215  }
216  if (Ftrial - state_->value > c1_*Qk || Ftrial > FtrialP) {
217  state_->iterateVec->set(*pxP);
218  strial = strialP;
219  ntrial = ntrialP;
220  Ftrial = FtrialP;
221  state_->searchSize = alphaP;
222  state_->stepVec->set(*state_->iterateVec);
223  state_->stepVec->axpy(-one,x);
224  accept = false;
225  }
226  }
227  else {
228  while ( Ftrial - state_->value > c1_*Qk && ls_nfval < maxit_ ) {
229  // Decrease search size
230  state_->searchSize *= rhodec_;
231  // Compute proximal gradient step with new search size
232  pgstep(*state_->iterateVec, *state_->stepVec, nobj, x, *dg, state_->searchSize, tol);
233  // Compute objective function values
234  sobj.update(*state_->iterateVec,UpdateType::Trial);
235  strial = sobj.value(*state_->iterateVec,tol);
236  nobj.update(*state_->iterateVec,UpdateType::Trial);
237  ntrial = nobj.value(*state_->iterateVec,tol);
238  Ftrial = strial + ntrial;
239  ls_nfval++;
240  // Compute decrease indicator
241  Qk = state_->gradientVec->apply(*state_->stepVec) + ntrial - state_->nvalue;
242  if (verbosity_ > 1) {
243  outStream << std::endl;
244  outStream << " Step size: " << state_->searchSize << std::endl;
245  outStream << " Trial smooth value: " << strial << std::endl;
246  outStream << " Trial nonsmooth value: " << ntrial << std::endl;
247  outStream << " Computed reduction: " << state_->value-Ftrial << std::endl;
248  outStream << " Dot product of gradient and step: " << Qk << std::endl;
249  outStream << " Sufficient decrease bound: " << -Qk*c1_ << std::endl;
250  outStream << " Number of function evaluations: " << ls_nfval << std::endl;
251  }
252  }
253  }
254  state_->nsval += ls_nfval;
255  state_->nnval += ls_nfval;
256 
257  // Compute norm of step
258  state_->snorm = state_->stepVec->norm();
259 
260  // Update iterate
261  state_->iter++;
262  x.set(*state_->iterateVec);
263 
264  // Compute new value and gradient
265  state_->svalue = strial;
266  state_->nvalue = ntrial;
267  state_->value = Ftrial;
268  if (accept) {
269  sobj.update(x,UpdateType::Accept,state_->iter);
270  nobj.update(x,UpdateType::Accept,state_->iter);
271  }
272  else {
273  sobj.update(x,UpdateType::Revert,state_->iter);
274  nobj.update(x,UpdateType::Revert,state_->iter);
275  }
276  sobj.gradient(*state_->gradientVec,x,tol);
277  state_->ngrad++;
278  dg->set(state_->gradientVec->dual());
279 
280  // Compute proximal gradient step with initial search size
281  searchSize = state_->searchSize;
282  if (!usePrevAlpha_ && !useAdapt_) searchSize = alpha0_;
283  pgstep(*state_->iterateVec, *state_->stepVec, nobj, x, *dg, searchSize, tol);
284  snorm = state_->stepVec->norm();
285  state_->gnorm = snorm / searchSize;
286 
287  // Update Output
288  if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
289  }
290  if (verbosity_ > 0) TypeP::Algorithm<Real>::writeExitStatus(outStream);
291 }
292 
293 template<typename Real>
294 void ProxGradientAlgorithm<Real>::writeHeader( std::ostream& os ) const {
295  std::ios_base::fmtflags osFlags(os.flags());
296  if (verbosity_ > 1) {
297  os << std::string(109,'-') << std::endl;
298  os << "Proximal gradient descent";
299  os << " status output definitions" << std::endl << std::endl;
300  os << " iter - Number of iterates (steps taken)" << std::endl;
301  os << " value - Objective function value" << std::endl;
302  os << " gnorm - Norm of the proximal gradient with parameter alpha" << std::endl;
303  os << " snorm - Norm of the step (update to optimization vector)" << std::endl;
304  os << " alpha - Line search step length" << std::endl;
305  os << " #sval - Cumulative number of times the smooth objective function was evaluated" << std::endl;
306  os << " #nval - Cumulative number of times the nonsmooth objective function was evaluated" << std::endl;
307  os << " #grad - Cumulative number of times the gradient was computed" << std::endl;
308  os << " #prox - Cumulative number of times the proximal operator was computed" << std::endl;
309  os << std::string(109,'-') << std::endl;
310  }
311 
312  os << " ";
313  os << std::setw(6) << std::left << "iter";
314  os << std::setw(15) << std::left << "value";
315  os << std::setw(15) << std::left << "gnorm";
316  os << std::setw(15) << std::left << "snorm";
317  os << std::setw(15) << std::left << "alpha";
318  os << std::setw(10) << std::left << "#sval";
319  os << std::setw(10) << std::left << "#nval";
320  os << std::setw(10) << std::left << "#grad";
321  os << std::setw(10) << std::left << "#nprox";
322  os << std::endl;
323  os.flags(osFlags);
324 }
325 
326 template<typename Real>
327 void ProxGradientAlgorithm<Real>::writeName( std::ostream& os ) const {
328  std::ios_base::fmtflags osFlags(os.flags());
329  os << std::endl << "Proximal Gradient Descent with Bidirectional Line Search (Type P)" << std::endl;
330  os.flags(osFlags);
331 }
332 
333 template<typename Real>
334 void ProxGradientAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
335  std::ios_base::fmtflags osFlags(os.flags());
336  os << std::scientific << std::setprecision(6);
337  if ( state_->iter == 0 ) writeName(os);
338  if ( write_header ) writeHeader(os);
339  if ( state_->iter == 0 ) {
340  os << " ";
341  os << std::setw(6) << std::left << state_->iter;
342  os << std::setw(15) << std::left << state_->value;
343  os << std::setw(15) << std::left << state_->gnorm;
344  os << std::setw(15) << std::left << "---";
345  os << std::setw(15) << std::left << "---";
346  os << std::setw(10) << std::left << state_->nsval;
347  os << std::setw(10) << std::left << state_->nnval;
348  os << std::setw(10) << std::left << state_->ngrad;
349  os << std::setw(10) << std::left << state_->nprox;
350  os << std::endl;
351  }
352  else {
353  os << " ";
354  os << std::setw(6) << std::left << state_->iter;
355  os << std::setw(15) << std::left << state_->value;
356  os << std::setw(15) << std::left << state_->gnorm;
357  os << std::setw(15) << std::left << state_->snorm;
358  os << std::setw(15) << std::left << state_->searchSize;
359  os << std::setw(10) << std::left << state_->nsval;
360  os << std::setw(10) << std::left << state_->nnval;
361  os << std::setw(10) << std::left << state_->ngrad;
362  os << std::setw(10) << std::left << state_->nprox;
363  os << std::endl;
364  }
365  os.flags(osFlags);
366 }
367 
368 } // namespace TypeP
369 } // namespace ROL
370 
371 #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::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::Objective::prox
virtual void prox(Vector< Real > &Pv, const Vector< Real > &v, Real t, Real &tol)
Definition: ROL_Objective.hpp:189
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::TypeP::ProxGradientAlgorithm::run
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &sobj, Objective< Real > &nobj, std::ostream &outStream=std::cout) override
Run algorithm on unconstrained problems (Type-U). This general interface supports the use of dual opt...
Definition: ROL_TypeP_ProxGradientAlgorithm_Def.hpp:129
ROL::TypeP::ProxGradientAlgorithm::initialize
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &sobj, Objective< Real > &nobj, Vector< Real > &px, Vector< Real > &dg, std::ostream &outStream=std::cout)
Definition: ROL_TypeP_ProxGradientAlgorithm_Def.hpp:76
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:74
ROL::StatusTest
Provides an interface to check status of optimization algorithms.
Definition: ROL_StatusTest.hpp:58
ROL::TypeP::ProxGradientAlgorithm::writeOutput
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
Definition: ROL_TypeP_ProxGradientAlgorithm_Def.hpp:334
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::TypeP::ProxGradientAlgorithm::writeHeader
void writeHeader(std::ostream &os) const override
Print iterate header.
Definition: ROL_TypeP_ProxGradientAlgorithm_Def.hpp:294
ROL::TypeP::ProxGradientAlgorithm::writeName
void writeName(std::ostream &os) const override
Print step name.
Definition: ROL_TypeP_ProxGradientAlgorithm_Def.hpp:327
ROL::TypeP::Algorithm::writeExitStatus
virtual void writeExitStatus(std::ostream &os) const
Definition: ROL_TypeP_Algorithm_Def.hpp:174
ROL::TypeP::Algorithm::initialize
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Definition: ROL_TypeP_Algorithm_Def.hpp:62
Generated by   doxygen 1.8.5