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