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