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