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 
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_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 
97 template<typename Real>
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 
113 template<typename Real>
114 void 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 
128 template<typename Real>
130  return xprim_->dot(x) + b_;
131 }
132 
133 template<typename Real>
134 void 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 
140 template<typename Real>
141 void 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
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
const Ptr< BoundConstraint< Real > > bnd_
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:80
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
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)
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
const Ptr< Constraint< Real > > con_
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Provides the interface to apply upper and lower bound constraints.
Real residual(const Vector< Real > &x) const
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
Defines the general constraint operator interface.