ROL
ROL_ColemanLiModel.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 #ifndef ROL_COLEMANLIMODEL_HPP
45 #define ROL_COLEMANLIMODEL_HPP
46 
47 #include "ROL_TrustRegionModel.hpp"
48 #include "ROL_BoundConstraint.hpp"
49 #include "ROL_Secant.hpp"
50 
59 namespace ROL {
60 
61 template<class Real>
62 class ColemanLiModel : public TrustRegionModel<Real> {
63 private:
64  Teuchos::RCP<BoundConstraint<Real> > bnd_; // Bound constraint
65  Teuchos::RCP<Secant<Real> > sec_; // Secant storage
66 
67  Teuchos::RCP<Vector<Real> > prim_, dual_, hv_; // Auxiliary storage
68  Teuchos::RCP<Vector<Real> > step_; // Step storage
69  Teuchos::RCP<Vector<Real> > cauchyStep_, cauchyScal_; // Cauchy point vectors
70  Teuchos::RCP<Vector<Real> > reflectStep_, reflectScal_; // Reflective step vectors
71  Teuchos::RCP<Vector<Real> > Dmat_; // sqrt(abs(v))
72  Teuchos::RCP<Vector<Real> > Cmat_; // diag(g) * dv/dx
73 
74  const bool useSecantPrecond_; // Use secant as preconditioner (unused)
75  const bool useSecantHessVec_; // Use secant as Hessian
76 
77  const Real TRradius_, stepBackMax_, stepBackScale_; // Primal transform parameters
78  const bool singleReflect_; // Use single reflection
79  Real sCs_, pred_; // Actual/predicted reduction
80 
81  Elementwise::Multiply<Real> mult_; // Elementwise multiply
82  Elementwise::Divide<Real> div_; // Elementwise division
83 
84  // Apply diagonal D matrix
85  void applyD( Vector<Real> &Dv, const Vector<Real> &v ) {
86  Dv.set(v);
87  Dv.applyBinary(div_,*Dmat_);
88  }
89 
90  // Apply inverse of diagonal D matrix
91  void applyInverseD( Vector<Real> &Dv, const Vector<Real> &v ) {
92  Dv.set(v);
93  Dv.applyBinary(mult_,*Dmat_);
94  }
95 
96  // Apply diagonal C matrix
97  void applyC( Vector<Real> &Cv, const Vector<Real> &v ) {
98  Cv.set(v);
99  Cv.applyBinary(mult_, *Cmat_);
100  }
101 
102  void constructC(void) {
103  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
104  const Teuchos::RCP<const Vector<Real> > l = bnd_->getLowerVectorRCP();
105  const Teuchos::RCP<const Vector<Real> > u = bnd_->getUpperVectorRCP();
106 
107  // Set Cmat_ to be the sign of the gradient
108  Cmat_->set(gc->dual());
109  Cmat_->applyUnary(Elementwise::Sign<Real>());
110  // If g < 0 and u = INF then set Cmat_ to zero
111  class NegGradInfU : public Elementwise::BinaryFunction<Real> {
112  public:
113  NegGradInfU(void) {}
114  Real apply(const Real &x, const Real &y) const {
115  const Real zero(0), one(1), INF(ROL_INF<Real>());
116  return (x < zero && y == INF) ? zero : one;
117  }
118  };
119  prim_->set(gc->dual());
120  prim_->applyBinary(NegGradInfU(), *u);
121  Cmat_->applyBinary(mult_, *prim_);
122  // If g >= 0 and l = -INF then set Cmat_ to zero
123  class PosGradNinfL : public Elementwise::BinaryFunction<Real> {
124  public:
125  PosGradNinfL(void) {}
126  Real apply(const Real &x, const Real &y) const {
127  const Real zero(0), one(1), NINF(ROL_NINF<Real>());
128  return (x >= zero && y == NINF) ? zero : one;
129  }
130  };
131  prim_->set(gc->dual());
132  prim_->applyBinary(PosGradNinfL(), *l);
133  Cmat_->applyBinary(mult_, *prim_);
134  // Pointwise multiply Cmat_ with the gradient
135  Cmat_->applyBinary(mult_, gc->dual());
136  }
137 
138  void constructInverseD(void) {
139  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
140  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
141  const Teuchos::RCP<const Vector<Real> > l = bnd_->getLowerVectorRCP();
142  const Teuchos::RCP<const Vector<Real> > u = bnd_->getUpperVectorRCP();
143  const Real zero(0), one(1), INF(ROL_INF<Real>()), NINF(ROL_NINF<Real>());
144  const int LESS_THAN = 0;
145  const int EQUAL_TO = 1;
146  const int GREATER_THAN = 2;
147 
148  Dmat_->zero();
149  // CASE (i)
150  // Mask for negative gradient (m1 is 1 if g is negative and 0 otherwise)
151  reflectStep_->applyBinary(Elementwise::ValueSet<Real>(zero, LESS_THAN),gc->dual());
152  // Mask for finite upper bounds (m2 is 1 if u is finite and 0 otherwise)
153  reflectScal_->applyBinary(Elementwise::ValueSet<Real>(INF, LESS_THAN),*u);
154  // Mask for g_i < 0 and u_i < inf
155  reflectScal_->applyBinary(mult_,*reflectStep_);
156  // prim_i = { u_i-x_i if g_i < 0 and u_i < inf
157  // { 0 otherwise
158  prim_->set(*u); prim_->axpy(-one,*xc);
159  prim_->applyBinary(mult_,*reflectScal_);
160  // Add to D
161  Dmat_->plus(*prim_);
162 
163  // CASE (iii)
164  // Mask for infinite upper bounds
165  reflectScal_->applyBinary(Elementwise::ValueSet<Real>(INF, EQUAL_TO),*u);
166  // Mask for g_i < 0 and u_i = inf
167  reflectScal_->applyBinary(mult_,*reflectStep_);
168  // prim_i = { -1 if g_i < 0 and u_i = inf
169  // { 0 otherwise
170  prim_->applyUnary(Elementwise::Fill<Real>(-one));
171  prim_->applyBinary(mult_,*reflectScal_);
172  // Add to D
173  Dmat_->plus(*prim_);
174 
175  // CASE (ii)
176  // m1 = 1-m1
177  reflectStep_->scale(-one);
178  reflectStep_->applyUnary(Elementwise::Shift<Real>(one));
179  // Mask for finite lower bounds
180  reflectScal_->applyBinary(Elementwise::ValueSet<Real>(NINF, GREATER_THAN),*l);
181  // Zero out elements of Jacobian with l_i=-inf
182  reflectScal_->applyBinary(mult_,*reflectStep_);
183  // prim_i = { x_i-l_i if g_i >= 0 and l_i > -inf
184  // { 0 otherwise
185  prim_->set(*xc); prim_->axpy(-one,*l);
186  prim_->applyBinary(mult_,*reflectScal_);
187  // Add to D
188  Dmat_->plus(*prim_);
189 
190  // CASE (iv)
191  // Mask for infinite lower bounds
192  reflectScal_->applyBinary(Elementwise::ValueSet<Real>(NINF, EQUAL_TO),*l);
193  // Mask for g_i>=0 and l_i=-inf
194  reflectScal_->applyBinary(mult_,*reflectStep_);
195  // prim_i = { 1 if g_i >= 0 and l_i = -inf
196  // { 0 otherwise
197  prim_->applyUnary(Elementwise::Fill<Real>(one));
198  prim_->applyBinary(mult_,*reflectScal_);
199  // Add to D
200  Dmat_->plus(*prim_);
201 
202  // d_i = { u_i-x_i if g_i < 0, u_i<inf
203  // { -1 if g_i < 0, u_i=inf
204  // { x_i-l_i if g_i >= 0, l_i>-inf
205  // { 1 if g_i >= 0, l_i=-inf
206  Dmat_->applyUnary(Elementwise::AbsoluteValue<Real>());
207  Dmat_->applyUnary(Elementwise::SquareRoot<Real>());
208  }
209 
210  // Build diagonal D and C matrices
212  const Vector<Real> &x, const Vector<Real> &g) {
213  bnd_ = Teuchos::rcpFromRef(bnd);
214 
215  prim_ = x.clone();
216  dual_ = g.clone();
217  hv_ = g.clone();
218  step_ = x.clone();
219  Dmat_ = x.clone();
220  Cmat_ = x.clone();
221 
222  cauchyStep_ = x.clone();
223  cauchyScal_ = x.clone();
224  reflectStep_ = x.clone();
225  reflectScal_ = x.clone();
226 
227  constructC();
229  }
230 
231  public:
232 
234  const Vector<Real> &x, const Vector<Real> &g)
235  : TrustRegionModel<Real>::TrustRegionModel(obj,x,g,false),
236  sec_(Teuchos::null), useSecantPrecond_(false), useSecantHessVec_(false),
237  TRradius_(1), stepBackMax_(0.9999), stepBackScale_(1),
238  singleReflect_(true), sCs_(0), pred_(0) {
239  initialize(obj,bnd,x,g);
240  }
241 
243  const Vector<Real> &x, const Vector<Real> &g,
244  const Real TRradius, const Real stepBackMax, const Real stepBackScale,
245  const bool singleReflect = true )
246  : TrustRegionModel<Real>::TrustRegionModel(obj,x,g,false),
247  sec_(Teuchos::null), useSecantPrecond_(false), useSecantHessVec_(false),
248  TRradius_(TRradius), stepBackMax_(stepBackMax), stepBackScale_(stepBackScale),
249  singleReflect_(singleReflect), sCs_(0), pred_(0) {
250  initialize(obj,bnd,x,g);
251  }
252 
254  const Vector<Real> &x, const Vector<Real> &g,
255  const Teuchos::RCP<Secant<Real> > &sec,
256  const bool useSecantPrecond, const bool useSecantHessVec,
257  const Real TRradius, const Real stepBackMax, const Real stepBackScale,
258  const bool singleReflect = true )
259  : TrustRegionModel<Real>::TrustRegionModel(obj,x,g,false),
260  sec_(sec), useSecantPrecond_(useSecantPrecond), useSecantHessVec_(useSecantHessVec),
261  TRradius_(TRradius), stepBackMax_(stepBackMax), stepBackScale_(stepBackScale),
262  singleReflect_(singleReflect), sCs_(0), pred_(0) {
263  initialize(obj,bnd,x,g);
264  }
265 
266  // Note that s is the \f$\hat{s}\f$ and \f$\psi\f$ is the $\hat\psi$ from the paper
267  Real value( const Vector<Real> &s, Real &tol ) {
268  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
269  // Apply Hessian to s
270  hessVec(*hv_, s, s, tol);
271  hv_->scale(static_cast<Real>(0.5));
272  // Form inv(D) * g
273  applyInverseD(*prim_, gc->dual());
274  // Add scaled gradient to Hessian in direction s
275  hv_->plus(prim_->dual());
276  return hv_->dot(s.dual());
277  }
278 
279  void gradient( Vector<Real> &g, const Vector<Real> &s, Real &tol ) {
280  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
281  hessVec(g, s, s, tol);
282  applyInverseD(*prim_, gc->dual());
283  g.plus(prim_->dual());
284  }
285 
286  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &s, Real &tol ) {
287  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
288  // Build B = inv(D) * Hessian * inv(D)
289  applyInverseD(*prim_, v);
290  if(useSecantHessVec_) {
291  sec_->applyB(*dual_, *prim_);
292  }
293  else {
294  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
295  TrustRegionModel<Real>::getObjective()->hessVec(*dual_, *prim_, *xc, tol);
296  }
297  applyInverseD(hv, *dual_);
298  // Build C = diag(g) J
299  applyC(*prim_, v);
300  hv.plus(prim_->dual());
301  }
302 
303  void dualTransform( Vector<Real> &tv, const Vector<Real> &v ) {
304  applyInverseD(tv, v);
305  }
306 
307  void primalTransform( Vector<Real> &tiv, const Vector<Real> &v ) {
308  Real tol = std::sqrt(ROL_EPSILON<Real>());
309 
310  /**************************************************************************/
311  /* PERFORM OPTIMAL SCALING OF TRUST REGION SUBPROBLEM SOLUTION */
312  /**************************************************************************/
313  applyInverseD(tiv, v);
314  // Get bounds on scalar variable
315  Real lowerBoundV(ROL_NINF<Real>()), upperBoundV(ROL_INF<Real>());
316  getScalarBounds(lowerBoundV, upperBoundV, tiv);
317  // Minimize one dimensional quadratic over bounds
318  Real tauV(1);
319  Real valueV = minimize1D(tauV, lowerBoundV, upperBoundV, v);
320 
321  /**************************************************************************/
322  /* COMPUTE CAUCHY POINT: STORED IN cauchyStep_ AND cauchyScal_ */
323  /**************************************************************************/
324  Real valueG = computeCauchyPoint();
325 
326  /**************************************************************************/
327  /* COMPUTE REFLECTIVE STEP: STORED IN reflectStep_ AND reflectScal_ */
328  /**************************************************************************/
329  if ( singleReflect_ ) {
331  }
332  else {
334  }
336  // Get bounds on scalar variable
337  Real lowerBoundR(ROL_NINF<Real>()), upperBoundR(ROL_INF<Real>());
338  getScalarBounds(lowerBoundR, upperBoundR, *reflectScal_);
339  // Minimize one dimensional quadratic over bounds
340  Real tauR(1);
341  Real valueR = minimize1D(tauR, lowerBoundR, upperBoundR, *reflectStep_);
342 
343  /**************************************************************************/
344  /* CHOOSE STEP THAT GIVES MOST PREDICTED REDUCTION */
345  /**************************************************************************/
346  Real VALUE(0);
347  bool useCauchyPoint = (valueG < valueV);
348  if (useCauchyPoint) {
349  VALUE = valueG;
350  tiv.set(*cauchyScal_);
351  // Store unscaled step
352  step_->set(*cauchyStep_);
353  }
354  else {
355  VALUE = valueV;
356  tiv.scale(tauV);
357  // Store unscaled step
358  step_->set(v);
359  step_->scale(tauV);
360  }
361  bool useReflectStep = (valueR < VALUE);
362  if (useReflectStep) {
363  VALUE = valueR;
364  tiv.set(*reflectScal_);
365  tiv.scale(tauR);
366  // Store unscaled step
367  step_->set(*reflectStep_);
368  step_->scale(tauR);
369  }
370 
371  /**************************************************************************/
372  /* ENSURE CHOSEN STEP IS STRICTLY FEASIBLE */
373  /**************************************************************************/
374  // Computed predicted reduction based on input step
375  if ( !isStrictlyFeasibleStep(tiv) ) {
376  Real snorm = step_->norm();
377  Real theta = std::max( stepBackMax_, static_cast<Real>(1) - stepBackScale_ * snorm);
378  tiv.scale(theta);
379  step_->scale(theta);
380  // Compute predicted reduction
381  pred_ = -value(*step_,tol);
382  }
383  else { // Theta is one
384  // Compute predicted reduction
385  pred_ = -VALUE;
386  }
387 
388  // Compute update for actual reduction
389  applyC(*prim_, *step_);
390  sCs_ = static_cast<Real>(-0.5) * prim_->dot(*step_);
391  }
392 
393  void updatePredictedReduction(Real &pred, const Vector<Real> &s) {
394  pred = pred_;
395  }
396 
397  void updateActualReduction(Real &ared, const Vector<Real> &s) {
398  ared += sCs_;
399  }
400 
401  const Teuchos::RCP<BoundConstraint<Real> > getBoundConstraint(void) const {
402  return bnd_;
403  }
404 
405 private:
406 
407  void getScalarBounds( Real &lowerBound, Real &upperBound, const Vector<Real> &p ) {
408  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
409  const Teuchos::RCP<const Vector<Real> > l = bnd_->getLowerVectorRCP();
410  const Teuchos::RCP<const Vector<Real> > u = bnd_->getUpperVectorRCP();
411  const Real one(1);
412  Real pnorm = p.norm();
413 
414  // Define elementwise functions
415  class PruneNegative : public Elementwise::BinaryFunction<Real> {
416  private:
417  const Real val_;
418  public:
419  PruneNegative( const Real val ) : val_(val) {}
420  Real apply(const Real &x, const Real &y) const {
421  return (y < static_cast<Real>(0)) ? x/y : val_;
422  }
423  };
424  class PrunePositive : public Elementwise::BinaryFunction<Real> {
425  private:
426  const Real val_;
427  public:
428  PrunePositive( const Real val ) : val_(val) {}
429  Real apply(const Real &x, const Real &y) const {
430  return (y > static_cast<Real>(0)) ? x/y : val_;
431  }
432  };
433 
434  // Max of (l-x)/p if p > 0
435  prim_->set(*l); prim_->axpy(-one,*xc);
436  prim_->applyBinary(PrunePositive(ROL_NINF<Real>()),p);
437  Real lowerBound1 = prim_->reduce(Elementwise::ReductionMax<Real>());
438  // Max of (u-x)/p if p < 0
439  prim_->set(*u); prim_->axpy(-one,*xc);
440  prim_->applyBinary(PruneNegative(ROL_NINF<Real>()),p);
441  Real lowerBound2 = prim_->reduce(Elementwise::ReductionMax<Real>());
442  // Lower bound
443  Real lowerBound3 = std::max(lowerBound1, lowerBound2);
444 
445  // Min of (u-x)/p if p > 0
446  prim_->set(*u); prim_->axpy(-one,*xc);
447  prim_->applyBinary(PrunePositive(ROL_INF<Real>()),p);
448  Real upperBound1 = prim_->reduce(Elementwise::ReductionMin<Real>());
449  // Max of (l-x)/p if p < 0
450  prim_->set(*l); prim_->axpy(-one,*xc);
451  prim_->applyBinary(PruneNegative(ROL_INF<Real>()),p);
452  Real upperBound2 = prim_->reduce(Elementwise::ReductionMin<Real>());
453  // Upper bound
454  Real upperBound3 = std::min(upperBound1, upperBound2);
455 
456  // Adjust for trust-region constraint
457  lowerBound = std::max(lowerBound3, -TRradius_/pnorm);
458  upperBound = std::min(upperBound3, TRradius_/pnorm);
459  }
460 
461  Real minimize1D(Real &tau, const Real lowerBound, const Real upperBound, const Vector<Real> &p) {
462  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
463  Real tol = std::sqrt(ROL_EPSILON<Real>());
464 
465  // Compute coefficients of one dimensional quadratic
466  hessVec(*hv_, p, p, tol);
467  Real c2 = static_cast<Real>(0.5) * hv_->dot(p.dual());
468  applyInverseD(*prim_, gc->dual());
469  Real c1 = prim_->dot(p);
470 
471  // Minimize one dimensional quadratic over bounds
472  Real lval = (c2 * lowerBound + c1) * lowerBound;
473  Real rval = (c2 * upperBound + c1) * upperBound;
474  tau = (lval < rval) ? lowerBound : upperBound;
475  if (c2 > static_cast<Real>(0)) {
476  Real uncMin = static_cast<Real>(-0.5) * c1/c2;
477  tau = (uncMin > lowerBound && uncMin < upperBound) ? uncMin : tau;
478  }
479 
480  // Return minimal function value
481  return (c2 * tau + c1) * tau;
482  }
483 
484  Real computeCauchyPoint(void) {
485  // Set step = -inv(D) g
486  const Teuchos::RCP<const Vector<Real> > gc = TrustRegionModel<Real>::getGradient();
487  applyInverseD(*cauchyStep_, gc->dual());
488  cauchyStep_->scale(static_cast<Real>(-1));
489 
490  // Scale Cauchy point
492 
493  // Scalar bounds
494  Real lowerBound(ROL_NINF<Real>()), upperBound(ROL_INF<Real>());
495  getScalarBounds(lowerBound, upperBound, *cauchyScal_);
496 
497  // Minimize 1D quadratic over bounds
498  Real tau(1), value(0);
499  value = minimize1D(tau, lowerBound, upperBound, *cauchyStep_);
500 
501  // Scale Cauchy point and return minimal function value
502  cauchyStep_->scale(tau);
503  cauchyScal_->scale(tau);
504  return value;
505  }
506 
508  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
509  Real alpha = computeAlpha(Dv);
510  Rv.set(v);
511 
512  class LowerBound : public Elementwise::BinaryFunction<Real> {
513  public:
514  Real apply( const Real &x, const Real &y ) const {
515  return (x == y) ? static_cast<Real>(-1) : static_cast<Real>(1);
516  }
517  };
518  prim_->set(*xc); prim_->axpy(alpha,Dv);
519  prim_->applyBinary(LowerBound(),*bnd_->getLowerVectorRCP());
520  Rv.applyBinary(mult_,*prim_);
521 
522  class UpperBound : public Elementwise::BinaryFunction<Real> {
523  public:
524  Real apply( const Real &x, const Real &y ) const {
525  return (x == y) ? static_cast<Real>(-1) : static_cast<Real>(1);
526  }
527  };
528  prim_->set(*xc); prim_->axpy(alpha,Dv);
529  prim_->applyBinary(UpperBound(),*bnd_->getUpperVectorRCP());
530  Rv.applyBinary(mult_,*prim_);
531  }
532 
534  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
535  Rv.set(v);
536 
537  class LowerBound : public Elementwise::BinaryFunction<Real> {
538  public:
539  Real apply( const Real &x, const Real &y ) const {
540  return (x < y) ? static_cast<Real>(-1) : static_cast<Real>(1);
541  }
542  };
543  prim_->set(*xc); prim_->plus(Dv);
544  prim_->applyBinary(LowerBound(),*bnd_->getLowerVectorRCP());
545  Rv.applyBinary(mult_,*prim_);
546 
547  class UpperBound : public Elementwise::BinaryFunction<Real> {
548  public:
549  Real apply( const Real &x, const Real &y ) const {
550  return (x > y) ? static_cast<Real>(-1) : static_cast<Real>(1);
551  }
552  };
553  prim_->set(*xc); prim_->plus(Dv);
554  prim_->applyBinary(UpperBound(),*bnd_->getUpperVectorRCP());
555  Rv.applyBinary(mult_,*prim_);
556  }
557 
558  Real computeAlpha( const Vector<Real> &d ) {
559  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
560  Teuchos::RCP<Vector<Real> > lx = xc->clone(), ux = xc->clone();
561  const Real one(1);
562 
563  // Define elementwise functions
564  class SafeDivide : public Elementwise::BinaryFunction<Real> {
565  private:
566  const Real val_;
567  public:
568  SafeDivide( const Real val ) : val_(val) {}
569  Real apply(const Real &x, const Real &y) const {
570  const Real zero(0);
571  return (y == zero) ? val_ : x/y;
572  }
573  };
574 
575  // (l - x) / d
576  lx->set(*bnd_->getLowerVectorRCP());
577  lx->axpy(-one, *xc);
578  lx->applyBinary(SafeDivide(ROL_INF<Real>()), d);
579 
580  // (u - x) / d
581  ux->set(*bnd_->getUpperVectorRCP());
582  ux->axpy(-one, *xc);
583  ux->applyBinary(SafeDivide(ROL_INF<Real>()), d);
584 
585  // max{ (l - x) / d, (u - x) / d }
586  lx->applyBinary(Elementwise::Max<Real>(),*ux);
587 
588  // min{ max{ (l - x) / d, (u - x) / d } }
589  return lx->reduce(Elementwise::ReductionMin<Real>());
590  }
591 
592  bool isStrictlyFeasibleStep( const Vector<Real> &d ) const {
593  const Teuchos::RCP<const Vector<Real> > xc = TrustRegionModel<Real>::getIterate();
594 
595  class Greater : public Elementwise::BinaryFunction<Real> {
596  public:
597  Greater() {}
598  Real apply(const Real &x, const Real &y) const {
599  return (x > y) ? static_cast<Real>(1) : static_cast<Real>(0);
600  }
601  };
602  prim_->set(*xc); prim_->plus(d);
603  prim_->applyBinary(Greater(),*bnd_->getLowerVectorRCP());
604  Real lowerFeasible = prim_->reduce(Elementwise::ReductionMin<Real>());
605 
606  class Lesser : public Elementwise::BinaryFunction<Real> {
607  public:
608  Lesser() {}
609  Real apply(const Real &x, const Real &y) const {
610  return (x < y) ? static_cast<Real>(1) : static_cast<Real>(0);
611  }
612  };
613  prim_->set(*xc); prim_->plus(d);
614  prim_->applyBinary(Lesser(),*bnd_->getUpperVectorRCP());
615  Real upperFeasible = prim_->reduce(Elementwise::ReductionMin<Real>());
616 
617  return (upperFeasible * lowerFeasible > 0);
618  }
619 
620 }; // class ColemanLiModel
621 
622 }
623 
624 #endif // ROL_COLEMANLIMODEL_HPP
Provides the interface to evaluate objective functions.
void initialize(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g)
virtual void scale(const Real alpha)=0
Compute where .
Teuchos::RCP< Vector< Real > > step_
virtual void plus(const Vector &x)=0
Compute , where .
void computeFullReflectiveStep(Vector< Real > &Rv, const Vector< Real > &v, const Vector< Real > &Dv)
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:222
void gradient(Vector< Real > &g, const Vector< Real > &s, Real &tol)
Compute gradient.
bool isStrictlyFeasibleStep(const Vector< Real > &d) const
void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
void updatePredictedReduction(Real &pred, const Vector< Real > &s)
void updateActualReduction(Real &ared, const Vector< Real > &s)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
void primalTransform(Vector< Real > &tiv, const Vector< Real > &v)
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
const Teuchos::RCP< BoundConstraint< Real > > getBoundConstraint(void) const
Teuchos::RCP< Vector< Real > > cauchyStep_
Provides the interface to evaluate trust-region model functions.
ColemanLiModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g, const Teuchos::RCP< Secant< Real > > &sec, const bool useSecantPrecond, const bool useSecantHessVec, const Real TRradius, const Real stepBackMax, const Real stepBackScale, const bool singleReflect=true)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Elementwise::Divide< Real > div_
Elementwise::Multiply< Real > mult_
Teuchos::RCP< Vector< Real > > Cmat_
Teuchos::RCP< BoundConstraint< Real > > bnd_
void applyD(Vector< Real > &Dv, const Vector< Real > &v)
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:213
Teuchos::RCP< Vector< Real > > reflectScal_
Teuchos::RCP< Secant< Real > > sec_
Real minimize1D(Real &tau, const Real lowerBound, const Real upperBound, const Vector< Real > &p)
Real value(const Vector< Real > &s, Real &tol)
Compute value.
virtual const Teuchos::RCP< const Vector< Real > > getGradient(void) const
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:70
Teuchos::RCP< Vector< Real > > reflectStep_
Teuchos::RCP< Vector< Real > > cauchyScal_
void applyInverseD(Vector< Real > &Dv, const Vector< Real > &v)
void applyC(Vector< Real > &Cv, const Vector< Real > &v)
ColemanLiModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g, const Real TRradius, const Real stepBackMax, const Real stepBackScale, const bool singleReflect=true)
Provides the interface to evaluate interior trust-region model functions from the Coleman-Li bound co...
void getScalarBounds(Real &lowerBound, Real &upperBound, const Vector< Real > &p)
Real computeAlpha(const Vector< Real > &d)
Provides the interface to apply upper and lower bound constraints.
Teuchos::RCP< Vector< Real > > hv_
Teuchos::RCP< Vector< Real > > prim_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
Teuchos::RCP< Vector< Real > > Dmat_
virtual Real norm() const =0
Returns where .
ColemanLiModel(Objective< Real > &obj, BoundConstraint< Real > &bnd, const Vector< Real > &x, const Vector< Real > &g)
Teuchos::RCP< Vector< Real > > dual_
void computeReflectiveStep(Vector< Real > &Rv, const Vector< Real > &v, const Vector< Real > &Dv)
virtual const Teuchos::RCP< Objective< Real > > getObjective(void) const
virtual const Teuchos::RCP< const Vector< Real > > getIterate(void) const