ROL
ROL_SemismoothNewtonProjection_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_SEMISMOOTHNEWTONPROJECTION_DEF_H
11 #define ROL_SEMISMOOTHNEWTONPROJECTION_DEF_H
12 
13 namespace ROL {
14 
15 template<typename Real>
17  const Vector<Real> &xdual,
18  const Ptr<BoundConstraint<Real>> &bnd,
19  const Ptr<Constraint<Real>> &con,
20  const Vector<Real> &mul,
21  const Vector<Real> &res)
22  : PolyhedralProjection<Real>(xprim,xdual,bnd,con,mul,res),
23  DEFAULT_atol_ (std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
24  DEFAULT_rtol_ (std::sqrt(ROL_EPSILON<Real>())),
25  DEFAULT_stol_ (std::sqrt(ROL_EPSILON<Real>())),
26  DEFAULT_decr_ (1e-4),
27  DEFAULT_factor_ (0.5),
28  DEFAULT_regscale_ (1e-4),
29  DEFAULT_errscale_ (1e-2),
30  DEFAULT_maxit_ (5000),
31  DEFAULT_lstype_ (0),
32  DEFAULT_verbosity_ (0),
33  DEFAULT_useproj_ (false),
34  atol_ (DEFAULT_atol_),
35  rtol_ (DEFAULT_rtol_),
36  stol_ (DEFAULT_stol_),
37  decr_ (DEFAULT_decr_),
38  factor_ (DEFAULT_factor_),
39  regscale_ (DEFAULT_regscale_),
40  errscale_ (DEFAULT_errscale_),
41  maxit_ (DEFAULT_maxit_),
42  lstype_ (DEFAULT_lstype_),
43  verbosity_ (DEFAULT_verbosity_),
44  useproj_ (DEFAULT_useproj_) {
45  dim_ = mul.dimension();
46  xnew_ = xprim.clone();
47  lnew_ = mul.clone();
48  dlam_ = mul.clone();
49 
50  ParameterList list;
51  list.sublist("General").sublist("Krylov").set("Type", "Conjugate Gradients");
52  list.sublist("General").sublist("Krylov").set("Absolute Tolerance", 1e-6);
53  list.sublist("General").sublist("Krylov").set("Relative Tolerance", 1e-4);
54  list.sublist("General").sublist("Krylov").set("Iteration Limit", dim_);
55  list.sublist("General").set("Inexact Hessian-Times-A-Vector", false);
56  krylov_ = KrylovFactory<Real>(list);
57 
59 }
60 
61 template<typename Real>
63  const Vector<Real> &xdual,
64  const Ptr<BoundConstraint<Real>> &bnd,
65  const Ptr<Constraint<Real>> &con,
66  const Vector<Real> &mul,
67  const Vector<Real> &res,
68  ParameterList &list)
69  : PolyhedralProjection<Real>(xprim,xdual,bnd,con,mul,res),
70  DEFAULT_atol_ (std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
71  DEFAULT_rtol_ (std::sqrt(ROL_EPSILON<Real>())),
72  DEFAULT_stol_ (std::sqrt(ROL_EPSILON<Real>())),
73  DEFAULT_decr_ (1e-4),
74  DEFAULT_factor_ (0.5),
75  DEFAULT_regscale_ (1e-4),
76  DEFAULT_errscale_ (1e-2),
77  DEFAULT_maxit_ (5000),
78  DEFAULT_lstype_ (0),
79  DEFAULT_verbosity_ (0),
80  DEFAULT_useproj_ (false),
81  atol_ (DEFAULT_atol_),
82  rtol_ (DEFAULT_rtol_),
83  stol_ (DEFAULT_stol_),
84  decr_ (DEFAULT_decr_),
85  factor_ (DEFAULT_factor_),
86  regscale_ (DEFAULT_regscale_),
87  errscale_ (DEFAULT_errscale_),
88  maxit_ (DEFAULT_maxit_),
89  lstype_ (DEFAULT_lstype_),
90  verbosity_ (DEFAULT_verbosity_),
91  useproj_ (DEFAULT_useproj_) {
92  dim_ = mul.dimension();
93  xnew_ = xprim.clone();
94  lnew_ = mul.clone();
95  dlam_ = mul.clone();
96  ParameterList &ppl = list.sublist("General").sublist("Polyhedral Projection");
97  atol_ = ppl.get("Absolute Tolerance", DEFAULT_atol_);
98  rtol_ = ppl.get("Relative Tolerance", DEFAULT_rtol_);
99  stol_ = ppl.sublist("Semismooth Newton").get("Step Tolerance", DEFAULT_stol_);
100  decr_ = ppl.sublist("Semismooth Newton").get("Sufficient Decrease Tolerance", DEFAULT_decr_);
101  factor_ = ppl.sublist("Semismooth Newton").get("Backtracking Rate", DEFAULT_factor_);
102  regscale_ = ppl.sublist("Semismooth Newton").get("Regularization Scale", DEFAULT_regscale_);
103  errscale_ = ppl.sublist("Semismooth Newton").get("Relative Error Scale", DEFAULT_errscale_);
104  maxit_ = ppl.get("Iteration Limit", DEFAULT_maxit_);
105  lstype_ = ppl.sublist("Semismooth Newton").get("Line Search Type", DEFAULT_lstype_);
106  verbosity_ = list.sublist("General").get("Output Level", DEFAULT_verbosity_);
107  useproj_ = ppl.sublist("Semismooth Newton").get("Project onto Separating Hyperplane", DEFAULT_useproj_);
108 
109  ParameterList klist;
110  klist.sublist("General").sublist("Krylov") = ppl.sublist("Semismooth Newton").sublist("Krylov");
111  klist.sublist("General").set("Inexact Hessian-Times-A-Vector", false);
112  krylov_ = KrylovFactory<Real>(klist);
113 
115 }
116 
117 template<typename Real>
119  if (con_ == nullPtr) {
120  bnd_->project(x);
121  }
122  else {
123  project_ssn(x, *mul_, *dlam_, stream);
124  }
125 }
126 
127 template<typename Real>
129  Real tol(std::sqrt(ROL_EPSILON<Real>()));
130  con_->update(y,UpdateType::Temp);
131  con_->value(r,y,tol);
132  return r.norm();
133 }
134 
135 template<typename Real>
137  const Vector<Real> &r,
138  const Vector<Real> &y,
139  const Real mu,
140  const Real rho,
141  int &iter,
142  int &flag) const {
143  Ptr<Precond> M = makePtr<Precond>(mu);
144  Ptr<Jacobian> J = makePtr<Jacobian>(con_,bnd_,makePtrFromRef(y),xdual_,xprim_,mu);
145  krylov_->run(s,*J,r,*M,iter,flag);
146 }
147 
148 template<typename Real>
150  const Vector<Real> &x,
151  const Vector<Real> &lam) const {
152  Real tol(std::sqrt(ROL_EPSILON<Real>()));
153  y.set(x);
154  con_->update(x,UpdateType::Temp);
155  con_->applyAdjointJacobian(*xdual_,lam,x,tol);
156  y.plus(xdual_->dual());
157  bnd_->project(y);
158 }
159 
160 template<typename Real>
162  Vector<Real> &lam,
163  Vector<Real> &dlam,
164  std::ostream &stream) const {
165  const Real zero(0), half(0.5), one(1);
166  // Compute initial residual
167  update_primal(*xnew_,x,lam);
168  Real rnorm = residual(*res_,*xnew_);
169  if (rnorm == zero) {
170  x.set(*xnew_);
171  return;
172  }
173  Real alpha(1), tmp(0), mu(0), rho(1), dd(0);
174  int iter(0), flag(0);
175  std::ios_base::fmtflags streamFlags(stream.flags());
176  if (verbosity_ > 2) {
177  stream << std::endl;
178  stream << std::scientific << std::setprecision(6);
179  stream << " Polyhedral Projection using Dual Semismooth Newton" << std::endl;
180  stream << " ";
181  stream << std::setw(6) << std::left << "iter";
182  stream << std::setw(15) << std::left << "rnorm";
183  stream << std::setw(15) << std::left << "alpha";
184  stream << std::setw(15) << std::left << "mu";
185  stream << std::setw(15) << std::left << "rho";
186  stream << std::setw(15) << std::left << "rtol";
187  stream << std::setw(8) << std::left << "kiter";
188  stream << std::setw(8) << std::left << "kflag";
189  stream << std::endl;
190  }
191  for (int cnt = 0; cnt < maxit_; ++cnt) {
192  // Compute Newton step
193  mu = regscale_*std::max(rnorm,std::sqrt(rnorm));
194  rho = std::min(half,errscale_*std::min(std::sqrt(rnorm),rnorm));
195  solve_newton_system(dlam,*res_,*xnew_,mu,rho,iter,flag);
196  lnew_->set(lam); lnew_->axpy(-alpha, dlam);
197  update_primal(*xnew_,x,*lnew_);
198  // Perform line search
199  if (lstype_ == 1) { // Usual line search for nonlinear equations
200  tmp = residual(*res_,*xnew_);
201  while ( tmp > (one-decr_*alpha)*rnorm && alpha > stol_ ) {
202  alpha *= factor_;
203  lnew_->set(lam); lnew_->axpy(-alpha, dlam);
204  update_primal(*xnew_,x,*lnew_);
205  tmp = residual(*res_,*xnew_);
206  }
207  rnorm = tmp;
208  }
209  else { // Default Solodov and Svaiter line search
210  rnorm = residual(*res_,*xnew_);
211  //tmp = dlam.dot(res_->dual());
212  tmp = dlam.apply(*res_);
213  dd = dlam.dot(dlam);
214  while ( tmp < decr_*(one-rho)*mu*dd && alpha > stol_ ) {
215  alpha *= factor_;
216  lnew_->set(lam); lnew_->axpy(-alpha, dlam);
217  update_primal(*xnew_,x,*lnew_);
218  rnorm = residual(*res_,*xnew_);
219  //tmp = dlam.dot(res_->dual());
220  tmp = dlam.apply(*res_);
221  }
222  }
223  // Update iterate
224  lam.set(*lnew_);
225  // Project onto separating hyperplane
226  if (useproj_) {
227  lam.axpy(-alpha*tmp/(rnorm*rnorm),res_->dual());
228  update_primal(*xnew_,x,lam);
229  rnorm = residual(*res_,*xnew_);
230  }
231  if (verbosity_ > 2) {
232  stream << " ";
233  stream << std::setw(6) << std::left << cnt;
234  stream << std::setw(15) << std::left << rnorm;
235  stream << std::setw(15) << std::left << alpha;
236  stream << std::setw(15) << std::left << mu;
237  stream << std::setw(15) << std::left << rho;
238  stream << std::setw(15) << std::left << ctol_;
239  stream << std::setw(8) << std::left << iter;
240  stream << std::setw(8) << std::left << flag;
241  stream << std::endl;
242  }
243  if (rnorm <= ctol_) break;
244  alpha = one;
245  }
246  if (verbosity_ > 2) {
247  stream << std::endl;
248  }
249  if (rnorm > ctol_) {
250  //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
251  stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
252  stream << rnorm << " rtol = " << ctol_ << std::endl;
253  }
254  x.set(*xnew_);
255  stream.flags(streamFlags);
256 }
257 
258 template<typename Real>
260  // Set tolerance
261  Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
262  if (bnd_->isLowerActivated()) resl = residual(*res_,*bnd_->getLowerBound());
263  if (bnd_->isUpperActivated()) resu = residual(*res_,*bnd_->getUpperBound());
264  Real res0 = std::max(resl,resu);
265  if (res0 < atol_) res0 = static_cast<Real>(1);
266  return std::min(atol_,rtol_*res0);
267 }
268 
269 } // namespace ROL
270 
271 #endif
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:204
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
SemismoothNewtonProjection(const Vector< Real > &xprim, const Vector< Real > &xdual, const Ptr< BoundConstraint< Real >> &bnd, const Ptr< Constraint< Real >> &con, const Vector< Real > &mul, const Vector< Real > &res)
void solve_newton_system(Vector< Real > &s, const Vector< Real > &r, const Vector< Real > &y, const Real mu, const Real rho, int &iter, int &flag) const
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides the interface to apply upper and lower bound constraints.
Real residual(Vector< Real > &r, const Vector< Real > &y) const
void project_ssn(Vector< Real > &x, Vector< Real > &lam, Vector< Real > &dlam, std::ostream &stream=std::cout) const
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:57
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Vector< Real > &lam) const
Defines the general constraint operator interface.