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>
51 DaiFletcherProjection<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  initialize(xprim,xdual,bnd,con,mul,res);
69 }
70 
71 template<typename Real>
72 DaiFletcherProjection<Real>::DaiFletcherProjection(const Vector<Real> &xprim,
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 
98 template<typename Real>
99 void DaiFletcherProjection<Real>::initialize(const Vector<Real> &xprim,
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 
134 template<typename Real>
135 void DaiFletcherProjection<Real>::project(Vector<Real> &x, std::ostream &stream) {
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 
151 template<typename Real>
152 Real DaiFletcherProjection<Real>::residual(const Vector<Real> &x) const {
153  return xprim_->dot(x) + b_;
154 }
155 
156 template<typename Real>
157 void DaiFletcherProjection<Real>::update_primal(Vector<Real> &y, const Vector<Real> &x, const Real lam) const {
158  y.set(x);
159  y.axpy(lam,*xprim_);
160  bnd_->project(y);
161 }
162 
163 template<typename Real>
164 void 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
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::DaiFletcherProjection::residual
Real residual(const Vector< Real > &x) const
Definition: ROL_DaiFletcherProjection_Def.hpp:152
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:145
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::project
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Definition: ROL_DaiFletcherProjection_Def.hpp:135
ROL::DaiFletcherProjection::initialize
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)
Definition: ROL_DaiFletcherProjection_Def.hpp:99
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:164
ROL::DaiFletcherProjection::update_primal
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Definition: ROL_DaiFletcherProjection_Def.hpp:157
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:73
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:86
Generated by   doxygen 1.8.5