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