ROL
ROL_BrentsProjection_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_BRENTSPROJECTION_DEF_H
11 #define ROL_BRENTSPROJECTION_DEF_H
12 
13 namespace ROL {
14 
15 template<typename Real>
16 BrentsProjection<Real>::BrentsProjection(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 BrentsProjection<Real>::BrentsProjection(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 BrentsProjection<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::BrentsProjection : 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 BrentsProjection<Real>::project(Vector<Real> &x, std::ostream &stream) {
101  if (con_ == nullPtr) {
102  bnd_->project(x);
103  }
104  else {
105  mul1_ = -residual(x)/cdot_;
106  //mul1_ = static_cast<Real>(0);
107  dlam1_ = static_cast<Real>(2);
108  //dlam1_ = static_cast<Real>(1)+std::abs(mul1_);
109  project_df(x, mul1_, dlam1_, stream);
110  mul_->setScalar(mul1_);
111  }
112 }
113 
114 template<typename Real>
115 Real BrentsProjection<Real>::residual(const Vector<Real> &x) const {
116  return xprim_->dot(x) + b_;
117 }
118 
119 template<typename Real>
120 void BrentsProjection<Real>::update_primal(Vector<Real> &y, const Vector<Real> &x, const Real lam) const {
121  y.set(x);
122  y.axpy(lam,*xprim_);
123  bnd_->project(y);
124 }
125 
126 template<typename Real>
127 void BrentsProjection<Real>::project_df(Vector<Real> &x, Real &lam, Real &dlam, std::ostream &stream) const {
128  const Real zero(0), one(1), c1(0.1);
129  Real lamLower(0), lamUpper(0), res(0), resLower(0), resUpper(0), s(0);
130  Real rtol = ctol_;
131  int cnt(0);
132  // Compute initial residual
133  update_primal(*xnew_,x,lam);
134  res = residual(*xnew_);
135  if (res == zero) {
136  x.set(*xnew_);
137  return;
138  }
139  std::ios_base::fmtflags streamFlags(stream.flags());
140  if (verbosity_ > 2) {
141  stream << std::scientific << std::setprecision(6);
142  stream << std::endl;
143  stream << " Polyhedral Projection using Brents' Algorithm" << std::endl;
144  stream << " Bracketing Phase" << std::endl;
145  }
146  // Bracketing phase
147  if ( res < zero ) {
148  lamLower = lam;
149  resLower = res;
150  lam += dlam;
151  update_primal(*xnew_,x,lam);
152  res = residual(*xnew_);
153  if (verbosity_ > 2) {
154  stream << " ";
155  stream << std::setw(6) << std::left << "iter";
156  stream << std::setw(15) << std::left << "lam";
157  stream << std::setw(15) << std::left << "res";
158  stream << std::setw(15) << std::left << "lower lam";
159  stream << std::setw(15) << std::left << "lower res";
160  stream << std::endl;
161  stream << " ";
162  stream << std::setw(6) << std::left << cnt;
163  stream << std::setw(15) << std::left << lam;
164  stream << std::setw(15) << std::left << res;
165  stream << std::setw(15) << std::left << lamLower;
166  stream << std::setw(15) << std::left << resLower;
167  stream << std::endl;
168  }
169  while ( res < zero && std::abs(res) > rtol && cnt < maxit_ ) {
170  s = std::max(resLower/res-one,c1);
171  dlam += dlam/s;
172  lamLower = lam;
173  resLower = res;
174  lam += dlam;
175  update_primal(*xnew_,x,lam);
176  res = residual(*xnew_);
177  cnt++;
178  if (verbosity_ > 2) {
179  stream << " ";
180  stream << std::setw(6) << std::left << cnt;
181  stream << std::setw(15) << std::left << lam;
182  stream << std::setw(15) << std::left << res;
183  stream << std::setw(15) << std::left << lamLower;
184  stream << std::setw(15) << std::left << resLower;
185  stream << std::endl;
186  }
187  }
188  lamUpper = lam;
189  resUpper = res;
190  }
191  else {
192  lamUpper = lam;
193  resUpper = res;
194  lam -= dlam;
195  update_primal(*xnew_,x,lam);
196  res = residual(*xnew_);
197  if (verbosity_ > 2) {
198  stream << " ";
199  stream << std::setw(6) << std::left << "iter";
200  stream << std::setw(15) << std::left << "lam";
201  stream << std::setw(15) << std::left << "res";
202  stream << std::setw(15) << std::left << "upper lam";
203  stream << std::setw(15) << std::left << "upper res";
204  stream << std::endl;
205  stream << " ";
206  stream << std::setw(6) << std::left << cnt;
207  stream << std::setw(15) << std::left << lam;
208  stream << std::setw(15) << std::left << res;
209  stream << std::setw(15) << std::left << lamUpper;
210  stream << std::setw(15) << std::left << resUpper;
211  stream << std::endl;
212  }
213  while ( res > zero && std::abs(res) > rtol && cnt < maxit_ ) {
214  s = std::max(resUpper/res-one,c1);
215  dlam += dlam/s;
216  lamUpper = lam;
217  resUpper = res;
218  lam -= dlam;
219  update_primal(*xnew_,x,lam);
220  res = residual(*xnew_);
221  cnt++;
222  if (verbosity_ > 2) {
223  stream << " ";
224  stream << std::setw(6) << std::left << cnt;
225  stream << std::setw(15) << std::left << lam;
226  stream << std::setw(15) << std::left << res;
227  stream << std::setw(15) << std::left << lamUpper;
228  stream << std::setw(15) << std::left << resUpper;
229  stream << std::endl;
230  }
231  }
232  lamLower = lam;
233  resLower = res;
234  }
235  if (verbosity_ > 2) {
236  stream << " Bracket: ";
237  stream << std::setw(15) << std::left << lamLower;
238  stream << std::setw(15) << std::left << lamUpper;
239  stream << std::endl;
240  }
241 
242  // Secant phase
243  //rtol = ctol_*std::max(one,std::min(std::abs(resLower),std::abs(resUpper)));
244  cnt = 0;
245  if (verbosity_ > 2) {
246  stream << std::endl;
247  stream << " Brents' Phase" << std::endl;
248  stream << " ";
249  stream << std::setw(6) << std::left << "iter";
250  stream << std::setw(15) << std::left << "rtol";
251  stream << std::setw(15) << std::left << "lam";
252  stream << std::setw(15) << std::left << "res";
253  stream << std::setw(15) << std::left << "lam low";
254  stream << std::setw(15) << std::left << "res low";
255  stream << std::setw(15) << std::left << "lam up";
256  stream << std::setw(15) << std::left << "res up";
257  stream << std::endl;
258  }
259  const Real half(0.5), two(2), three(3);
260  const Real eps(ROL_EPSILON<Real>()), tol0(rtol); // tol0(1e1*eps);
261  Real d1(1), d2(1), tol(1);
262  Real p(0), q(0), r(0), m(0);
263  lam = lamUpper; res = resUpper;
264  update_primal(*xnew_,x,lamUpper);
265  for (cnt = 0; cnt < maxit_; cnt++) {
266  if ((resUpper > zero && res > zero) || (resUpper <= zero && res <= zero)) {
267  lam = lamLower; res = resLower;
268  d1 = lamUpper-lamLower; d2 = d1;
269  }
270  if (std::abs(res) < std::abs(resUpper)) {
271  lamLower = lamUpper; lamUpper = lam; lam = lamLower;
272  resLower = resUpper; resUpper = res; res = resLower;
273  }
274  tol = two*eps*std::abs(lamUpper) + half*tol0;
275  m = half*(lam - lamUpper);
276  if (std::abs(m) <= tol || std::abs(resUpper) <= rtol) break;
277  if (std::abs(d2) < tol || std::abs(resLower) <= std::abs(resUpper)) {
278  d1 = m; d2 = d1;
279  }
280  else {
281  s = resUpper/resLower;
282  if (lamLower == lam) {
283  p = two*m*s;
284  q = one-s;
285  }
286  else {
287  q = resLower/res;
288  r = resUpper/res;
289  p = s*(two*m*q*(q-r)-(lamUpper-lamLower)*(r-one));
290  q = (q-one)*(r-one)*(s-one);
291  }
292  if (p > zero) q = -q;
293  else p = -p;
294  if (two*p < three*m*q-std::abs(tol*q) && p < std::abs(half*d2*q)) {
295  d2 = d1; d1 = p/q;
296  }
297  else {
298  d1 = m; d2 = d1;
299  }
300  }
301  lamLower = lamUpper; resLower = resUpper;
302  if (std::abs(d1) > tol) lamUpper += d1;
303  else if (m > zero) lamUpper += tol;
304  else lamUpper -= tol;
305  update_primal(*xnew_,x,lamUpper);
306  resUpper = residual(*xnew_);
307 
308  if (verbosity_ > 2) {
309  stream << " ";
310  stream << std::setw(6) << std::left << cnt;
311  stream << std::setw(15) << std::left << rtol;
312  stream << std::setw(15) << std::left << lam;
313  stream << std::setw(15) << std::left << res;
314  stream << std::setw(15) << std::left << lamLower;
315  stream << std::setw(15) << std::left << resLower;
316  stream << std::setw(15) << std::left << lamUpper;
317  stream << std::setw(15) << std::left << resUpper;
318  stream << std::endl;
319  }
320  }
321  if (verbosity_ > 2) {
322  if (cnt < maxit_) {
323  stream << " ";
324  stream << std::setw(6) << std::left << cnt;
325  stream << std::setw(15) << std::left << rtol;
326  stream << std::setw(15) << std::left << lam;
327  stream << std::setw(15) << std::left << res;
328  stream << std::setw(15) << std::left << lamLower;
329  stream << std::setw(15) << std::left << resLower;
330  stream << std::setw(15) << std::left << lamUpper;
331  stream << std::setw(15) << std::left << resUpper;
332  stream << std::endl;
333  }
334  stream << std::endl;
335  }
336  // Return projection
337  res = resUpper;
338  x.set(*xnew_);
339  if (std::abs(res) > rtol ) {
340  //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
341  stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
342  stream << std::abs(res) << " rtol = " << rtol << std::endl;
343  }
344  stream.flags(streamFlags);
345 }
346 
347 } // namespace ROL
348 
349 #endif
ROL::BrentsProjection::residual
Real residual(const Vector< Real > &x) const
Definition: ROL_BrentsProjection_Def.hpp:115
ROL::BrentsProjection::project_df
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Definition: ROL_BrentsProjection_Def.hpp:127
ROL::BrentsProjection::BrentsProjection
BrentsProjection(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_BrentsProjection_Def.hpp:16
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::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::BrentsProjection::project
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Definition: ROL_BrentsProjection_Def.hpp:100
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:39
ROL::BrentsProjection::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_BrentsProjection_Def.hpp:64
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
ROL::BrentsProjection::update_primal
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Definition: ROL_BrentsProjection_Def.hpp:120
Generated by   doxygen 1.8.5