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