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