ROL
ROL_DaiFletcherProjection_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_DAIFLETCHERPROJECTION_DEF_H
46#define ROL_DAIFLETCHERPROJECTION_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_ltol_ (ROL_EPSILON<Real>()),
61 DEFAULT_maxit_ (5000),
62 DEFAULT_verbosity_ (0),
63 atol_ (DEFAULT_atol_),
64 rtol_ (DEFAULT_rtol_),
65 ltol_ (DEFAULT_ltol_),
66 maxit_ (DEFAULT_maxit_),
67 verbosity_ (DEFAULT_verbosity_) {
68 initialize(xprim,xdual,bnd,con,mul,res);
69}
70
71template<typename Real>
73 const Vector<Real> &xdual,
74 const Ptr<BoundConstraint<Real>> &bnd,
75 const Ptr<Constraint<Real>> &con,
76 const Vector<Real> &mul,
77 const Vector<Real> &res,
78 ParameterList &list)
79 : PolyhedralProjection<Real>(xprim,xdual,bnd,con,mul,res),
80 DEFAULT_atol_ (std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
81 DEFAULT_rtol_ (std::sqrt(ROL_EPSILON<Real>())),
82 DEFAULT_ltol_ (ROL_EPSILON<Real>()),
83 DEFAULT_maxit_ (5000),
84 DEFAULT_verbosity_ (0),
85 atol_ (DEFAULT_atol_),
86 rtol_ (DEFAULT_rtol_),
87 ltol_ (DEFAULT_ltol_),
88 maxit_ (DEFAULT_maxit_),
89 verbosity_ (DEFAULT_verbosity_) {
90 atol_ = list.sublist("General").sublist("Polyhedral Projection").get("Absolute Tolerance", DEFAULT_atol_);
91 rtol_ = list.sublist("General").sublist("Polyhedral Projection").get("Relative Tolerance", DEFAULT_rtol_);
92 ltol_ = list.sublist("General").sublist("Polyhedral Projection").get("Multiplier Tolerance", DEFAULT_ltol_);
93 maxit_ = list.sublist("General").sublist("Polyhedral Projection").get("Iteration Limit", DEFAULT_maxit_);
94 verbosity_ = list.sublist("General").get("Output Level", DEFAULT_verbosity_);
95 initialize(xprim,xdual,bnd,con,mul,res);
96}
97
98template<typename Real>
100 const Vector<Real> &xdual,
101 const Ptr<BoundConstraint<Real>> &bnd,
102 const Ptr<Constraint<Real>> &con,
103 const Vector<Real> &mul,
104 const Vector<Real> &res) {
105 dim_ = mul.dimension();
106 ROL_TEST_FOR_EXCEPTION(dim_!=1,std::logic_error,
107 ">>> ROL::DaiFletcherProjection : The range of the linear constraint must be one dimensional!");
108 xnew_ = xprim.clone();
109 Px_ = xprim.clone();
110 mul1_ = static_cast<Real>(0);
111 dlam1_ = static_cast<Real>(2);
112 // con.value(x) = xprim_->dot(x) + b_
113 Real tol(std::sqrt(ROL_EPSILON<Real>()));
114 xprim_->zero();
115 con_->update(*xprim_,UpdateType::Temp);
116 con_->value(*res_,*xprim_,tol);
117 b_ = res_->dot(*res_->basis(0));
118 mul_->setScalar(static_cast<Real>(1));
119 con_->applyAdjointJacobian(*xdual_,*mul_,xprim,tol);
120 xprim_->set(xdual_->dual());
121 cdot_ = xprim_->dot(*xprim_);
122 // Set tolerance
123 //xnew_->zero();
124 //bnd_->project(*xnew_);
125 //Real res0 = std::abs(residual(*xnew_));
126 Real resl = ROL_INF<Real>(), resu = ROL_INF<Real>();
127 if (bnd_->isLowerActivated()) resl = residual(*bnd_->getLowerBound());
128 if (bnd_->isUpperActivated()) resu = residual(*bnd_->getUpperBound());
129 Real res0 = std::max(resl,resu);
130 if (res0 < atol_) res0 = static_cast<Real>(1);
131 ctol_ = std::min(atol_,rtol_*res0);
132}
133
134template<typename Real>
136 if (con_ == nullPtr) {
137 bnd_->project(x);
138 }
139 else {
140 Px_->set(x); bnd_->project(*Px_);
141 mul1_ = -residual(*Px_)/cdot_;
142 //mul1_ = -residual(x)/cdot_;
143 //mul1_ = static_cast<Real>(0);
144 dlam1_ = static_cast<Real>(2);
145 //dlam1_ = static_cast<Real>(1)+std::abs(mul1_);
146 project_df(x, mul1_, dlam1_, stream);
147 mul_->setScalar(mul1_);
148 }
149}
150
151template<typename Real>
153 return xprim_->dot(x) + b_;
154}
155
156template<typename Real>
158 y.set(x);
159 y.axpy(lam,*xprim_);
160 bnd_->project(y);
161}
162
163template<typename Real>
164void DaiFletcherProjection<Real>::project_df(Vector<Real> &x, Real &lam, Real &dlam, std::ostream &stream) const {
165 const Real zero(0), one(1), two(2), c1(0.1), c2(0.75), c3(0.25);
166 Real lamLower(0), lamUpper(0), lamNew(0), res(0), resLower(0), resUpper(0), s(0);
167 Real rtol = ctol_;
168 int cnt(0);
169 // Compute initial residual
170 update_primal(*xnew_,x,lam);
171 res = residual(*xnew_);
172 if (res == zero) {
173 x.set(*xnew_);
174 return;
175 }
176 std::ios_base::fmtflags streamFlags(stream.flags());
177 if (verbosity_ > 2) {
178 stream << std::scientific << std::setprecision(6);
179 stream << std::endl;
180 stream << " Polyhedral Projection using the Dai-Fletcher Algorithm" << std::endl;
181 stream << " Bracketing Phase" << std::endl;
182 }
183 // Bracketing phase
184 if ( res < zero ) {
185 lamLower = lam;
186 resLower = res;
187 lam += dlam;
188 update_primal(*xnew_,x,lam);
189 res = residual(*xnew_);
190 if (verbosity_ > 2) {
191 stream << " ";
192 stream << std::setw(6) << std::left << "iter";
193 stream << std::setw(15) << std::left << "lam";
194 stream << std::setw(15) << std::left << "res";
195 stream << std::setw(15) << std::left << "lower lam";
196 stream << std::setw(15) << std::left << "lower res";
197 stream << std::endl;
198 stream << " ";
199 stream << std::setw(6) << std::left << cnt;
200 stream << std::setw(15) << std::left << lam;
201 stream << std::setw(15) << std::left << res;
202 stream << std::setw(15) << std::left << lamLower;
203 stream << std::setw(15) << std::left << resLower;
204 stream << std::endl;
205 }
206 while ( res < zero && std::abs(res) > rtol && cnt < maxit_ ) {
207 s = std::max(resLower/res-one,c1);
208 dlam += dlam/s;
209 lamLower = lam;
210 resLower = res;
211 lam += dlam;
212 update_primal(*xnew_,x,lam);
213 res = residual(*xnew_);
214 cnt++;
215 if (verbosity_ > 2) {
216 stream << " ";
217 stream << std::setw(6) << std::left << cnt;
218 stream << std::setw(15) << std::left << lam;
219 stream << std::setw(15) << std::left << res;
220 stream << std::setw(15) << std::left << lamLower;
221 stream << std::setw(15) << std::left << resLower;
222 stream << std::endl;
223 }
224 }
225 lamUpper = lam;
226 resUpper = res;
227 }
228 else {
229 lamUpper = lam;
230 resUpper = res;
231 lam -= dlam;
232 update_primal(*xnew_,x,lam);
233 res = residual(*xnew_);
234 if (verbosity_ > 2) {
235 stream << " ";
236 stream << std::setw(6) << std::left << "iter";
237 stream << std::setw(15) << std::left << "lam";
238 stream << std::setw(15) << std::left << "res";
239 stream << std::setw(15) << std::left << "upper lam";
240 stream << std::setw(15) << std::left << "upper res";
241 stream << std::endl;
242 stream << " ";
243 stream << std::setw(6) << std::left << cnt;
244 stream << std::setw(15) << std::left << lam;
245 stream << std::setw(15) << std::left << res;
246 stream << std::setw(15) << std::left << lamUpper;
247 stream << std::setw(15) << std::left << resUpper;
248 stream << std::endl;
249 }
250 while ( res > zero && std::abs(res) > rtol && cnt < maxit_ ) {
251 s = std::max(resUpper/res-one,c1);
252 dlam += dlam/s;
253 lamUpper = lam;
254 resUpper = res;
255 lam -= dlam;
256 update_primal(*xnew_,x,lam);
257 res = residual(*xnew_);
258 cnt++;
259 if (verbosity_ > 2) {
260 stream << " ";
261 stream << std::setw(6) << std::left << cnt;
262 stream << std::setw(15) << std::left << lam;
263 stream << std::setw(15) << std::left << res;
264 stream << std::setw(15) << std::left << lamUpper;
265 stream << std::setw(15) << std::left << resUpper;
266 stream << std::endl;
267 }
268 }
269 lamLower = lam;
270 resLower = res;
271 }
272 if (verbosity_ > 2) {
273 stream << " Bracket: ";
274 stream << std::setw(15) << std::left << lamLower;
275 stream << std::setw(15) << std::left << lamUpper;
276 stream << std::endl;
277 }
278
279 // Secant phase
280 rtol = ctol_*std::max(one,std::min(std::abs(resLower),std::abs(resUpper)));
281 //s = one - resLower / resUpper;
282 //dlam = (lamUpper - lamLower) / s;
283 //lam = lamUpper - dlam;
284 s = (resUpper - resLower) / resUpper;
285 lam = (resUpper * lamLower - resLower * lamUpper) / (resUpper - resLower);
286 dlam = lamUpper - lam;
287 update_primal(*xnew_,x,lam);
288 res = residual(*xnew_);
289 cnt = 0;
290 if (verbosity_ > 2) {
291 stream << std::endl;
292 stream << " Secant Phase" << std::endl;
293 stream << " ";
294 stream << std::setw(6) << std::left << "iter";
295 stream << std::setw(15) << std::left << "lam";
296 stream << std::setw(15) << std::left << "res";
297 stream << std::setw(15) << std::left << "stepsize";
298 stream << std::setw(15) << std::left << "rtol";
299 stream << std::setw(15) << std::left << "lbnd";
300 stream << std::setw(15) << std::left << "lres";
301 stream << std::setw(15) << std::left << "ubnd";
302 stream << std::setw(15) << std::left << "ures";
303 stream << std::endl;
304 stream << " ";
305 stream << std::setw(6) << std::left << cnt;
306 stream << std::setw(15) << std::left << lam;
307 stream << std::setw(15) << std::left << res;
308 stream << std::setw(15) << std::left << dlam;
309 stream << std::setw(15) << std::left << rtol;
310 stream << std::setw(15) << std::left << lamLower;
311 stream << std::setw(15) << std::left << resLower;
312 stream << std::setw(15) << std::left << lamUpper;
313 stream << std::setw(15) << std::left << resUpper;
314 stream << std::endl;
315 }
316 for (cnt = 1; cnt < maxit_; cnt++) {
317 // Exit if residual or bracket length are sufficiently small
318 if ( std::abs(res) <= rtol ||
319 std::abs(lamUpper-lamLower) < ltol_*std::max(std::abs(lamUpper),std::abs(lamLower)) ) {
320 break;
321 }
322
323 if ( res > zero ) {
324 if ( s <= two ) {
325 lamUpper = lam;
326 resUpper = res;
327 //s = one - resLower / resUpper;
328 //dlam = (lamUpper - lamLower) / s;
329 //lam = lamUpper - dlam;
330 s = (resUpper - resLower) / resUpper;
331 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
332 dlam = lamUpper - lam;
333 }
334 else {
335 //s = std::max(resUpper / res - one, c1);
336 //dlam = (lamUpper - lam) / s;
337 //lamNew = std::max(lam - dlam, c2*lamLower + c3*lam);
338 if (resUpper <= (c1+one)*res) {
339 dlam = (lamUpper - lam) / c1;
340 lamNew = std::max(lam - dlam, c2*lamLower + c3*lam);
341 }
342 else {
343 lamNew = std::max((lam * resUpper - lamUpper * res) / (resUpper - res),
344 c2*lamLower + c3*lam);
345 dlam = lam - lamNew;
346 }
347 lamUpper = lam;
348 resUpper = res;
349 lam = lamNew;
350 s = (lamUpper - lamLower) / (lamUpper - lam);
351 }
352 }
353 else {
354 if ( s >= two ) {
355 lamLower = lam;
356 resLower = res;
357 //s = one - resLower / resUpper;
358 //dlam = (lamUpper - lamLower) / s;
359 //lam = lamUpper - dlam;
360 s = (resUpper - resLower) / resUpper;
361 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
362 dlam = lamUpper - lam;
363 }
364 else {
365 //s = std::max(resLower / res - one, c1);
366 //dlam = (lam + lamLower) / s;
367 //lamNew = std::min(lam + dlam, c2*lamUpper + c3*lam);
368 if (resLower >= (c1+one)*res) {
369 dlam = (lam - lamLower) / c1;
370 lamNew = std::max(lam + dlam, c2*lamUpper + c3*lam);
371 }
372 else {
373 lamNew = std::max((lamLower * res - lam * resLower) / (res - resLower),
374 c2*lamUpper + c3*lam);
375 dlam = lamNew - lamLower;
376 }
377 lamLower = lam;
378 resLower = res;
379 lam = lamNew;
380 s = (lamUpper - lamLower) / (lamUpper - lam);
381 }
382 }
383 update_primal(*xnew_,x,lam);
384 res = residual(*xnew_);
385
386 if (verbosity_ > 2) {
387 stream << " ";
388 stream << std::setw(6) << std::left << cnt;
389 stream << std::setw(15) << std::left << lam;
390 stream << std::setw(15) << std::left << res;
391 stream << std::setw(15) << std::left << dlam;
392 stream << std::setw(15) << std::left << rtol;
393 stream << std::setw(15) << std::left << lamLower;
394 stream << std::setw(15) << std::left << resLower;
395 stream << std::setw(15) << std::left << lamUpper;
396 stream << std::setw(15) << std::left << resUpper;
397 stream << std::endl;
398 }
399 }
400 if (verbosity_ > 2) {
401 stream << std::endl;
402 }
403 // Return projection
404 x.set(*xnew_);
405 if (std::abs(res) > rtol ) {
406 //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
407 stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
408 stream << std::abs(res) << " rtol = " << rtol << std::endl;
409 }
410 stream.flags(streamFlags);
411}
412
413} // namespace ROL
414
415#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 project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
DaiFletcherProjection(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)
Real residual(const Vector< Real > &x) const
void initialize(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 update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91