ROL
ROL_TrustRegionUtilities.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_TRUSTREGIONUTILITIES_U_H
45 #define ROL_TRUSTREGIONUTILITIES_U_H
46 
48 
49 namespace ROL {
50 namespace TRUtils {
51 
62 enum ETRFlag {
63  SUCCESS = 0,
70 };
71 
72 inline std::string ETRFlagToString(ETRFlag trf) {
73  std::string retString;
74  switch(trf) {
75  case SUCCESS:
76  retString = "Both actual and predicted reductions are positive (success)";
77  break;
78  case POSPREDNEG:
79  retString = "Actual reduction is positive and predicted reduction is negative (impossible)";
80  break;
81  case NPOSPREDPOS:
82  retString = "Actual reduction is nonpositive and predicted reduction is positive";
83  break;
84  case NPOSPREDNEG:
85  retString = "Actual reduction is nonpositive and predicted reduction is negative (impossible)";
86  break;
87  case TRNAN:
88  retString = "Actual and/or predicted reduction is a NaN";
89  break;
90  case QMINSUFDEC:
91  retString = "Subproblem solution did not produce sufficient decrease";
92  break;
93  default:
94  retString = "INVALID ETRFlag";
95  }
96  return retString;
97 }
98 
99 template<typename Real>
100 inline Real initialRadius(int &nfval,
101  const Vector<Real> &x,
102  const Vector<Real> &g,
103  Vector<Real> &Bg,
104  const Real fx,
105  const Real gnorm,
106  Objective<Real> &obj,
108  const Real delMax,
109  std::ostream &outStream,
110  const bool print = false) {
111  const Real zero(0), half(0.5), one(1), two(2), three(3), six(6);
112  const Real eps(ROL_EPSILON<Real>());
113  Real del(ROL_INF<Real>());
114  Ptr<Vector<Real>> xcp = x.clone();
115  model.setData(obj,x,g);
116  Real htol = std::sqrt(eps);
117  model.hessVec(Bg,g.dual(),x,htol);
118  Real gBg = Bg.dot(g);
119  Real alpha = one;
120  if ( gBg > eps ) {
121  alpha = gnorm*gnorm/gBg;
122  }
123  // Evaluate the objective function at the Cauchy point
124  xcp->set(g.dual());
125  xcp->scale(-alpha);
126  //Real gs = xcp->dot(g.dual());
127  Real gs = xcp->apply(g);
128  xcp->plus(x);
129  obj.update(*xcp,UpdateType::Temp);
130  Real ftol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>();
131  Real fnew = obj.value(*xcp,ftol); // MUST DO SOMETHING HERE WITH FTOL
132  nfval++;
133  // Perform cubic interpolation to determine initial trust region radius
134  Real a = fnew - fx - gs - half*alpha*alpha*gBg;
135  if ( std::abs(a) < eps ) {
136  // a = 0 implies the objective is quadratic in the negative gradient direction
137  del = std::min(alpha*gnorm,delMax);
138  }
139  else {
140  Real b = half*alpha*alpha*gBg;
141  Real c = gs;
142  if ( b*b-three*a*c > eps ) {
143  // There is at least one critical point
144  Real t1 = (-b-std::sqrt(b*b-three*a*c))/(three*a);
145  Real t2 = (-b+std::sqrt(b*b-three*a*c))/(three*a);
146  if ( six*a*t1 + two*b > zero ) {
147  // t1 is the minimizer
148  del = std::min(t1*alpha*gnorm,delMax);
149  }
150  else {
151  // t2 is the minimizer
152  del = std::min(t2*alpha*gnorm,delMax);
153  }
154  }
155  else {
156  del = std::min(alpha*gnorm,delMax);
157  }
158  }
159  if (del <= eps*gnorm) {
160  del = one;
161  }
162  obj.update(x,UpdateType::Revert);
163  if ( print ) {
164  outStream << " In TrustRegionUtilities::initialRadius" << std::endl;
165  outStream << " Initial radius: " << del << std::endl;
166  }
167  return del;
168 }
169 
170 template<typename Real>
171 inline void analyzeRatio(Real &rho,
172  ETRFlag &flag,
173  const Real fold,
174  const Real ftrial,
175  const Real pRed,
176  const Real epsi,
177  std::ostream &outStream = std::cout,
178  const bool print = false) {
179  const Real zero(0), one(1);
180  Real eps = epsi*std::max(one,fold);
181  Real aRed = fold - ftrial;
182  Real aRed_safe = aRed + eps, pRed_safe = pRed + eps;
183  if (((std::abs(aRed_safe) < epsi) && (std::abs(pRed_safe) < epsi)) || aRed == pRed) {
184  rho = one;
185  flag = SUCCESS;
186  }
187  else if ( std::isnan(aRed_safe) || std::isnan(pRed_safe) ) {
188  rho = -one;
189  flag = TRNAN;
190  }
191  else {
192  rho = aRed_safe/pRed_safe;
193  if (pRed_safe < zero && aRed_safe > zero) {
194  flag = POSPREDNEG;
195  }
196  else if (aRed_safe <= zero && pRed_safe > zero) {
197  flag = NPOSPREDPOS;
198  }
199  else if (aRed_safe <= zero && pRed_safe < zero) {
200  flag = NPOSPREDNEG;
201  }
202  else {
203  flag = SUCCESS;
204  }
205  }
206  if ( print ) {
207  outStream << " In TrustRegionUtilities::analyzeRatio" << std::endl;
208  outStream << " Current objective function value: " << fold << std::endl;
209  outStream << " New objective function value: " << ftrial << std::endl;
210  outStream << " Actual reduction: " << aRed << std::endl;
211  outStream << " Predicted reduction: " << pRed << std::endl;
212  outStream << " Safeguard: " << epsi << std::endl;
213  outStream << " Actual reduction with safeguard: " << aRed_safe << std::endl;
214  outStream << " Predicted reduction with safeguard: " << pRed_safe << std::endl;
215  outStream << " Ratio of actual and predicted reduction: " << rho << std::endl;
216  outStream << " Trust-region flag: " << flag << std::endl;
217  outStream << std::endl;
218  }
219 }
220 
221 template<typename Real>
222 inline Real interpolateRadius(const Vector<Real> &g,
223  const Vector<Real> &s,
224  const Real snorm,
225  const Real pRed,
226  const Real fold,
227  const Real ftrial,
228  const Real del,
229  const Real gamma0,
230  const Real gamma1,
231  const Real eta2,
232  std::ostream &outStream = std::cout,
233  const bool print = false) {
234  const Real one(1);
235  //Real gs = g.dot(s.dual());
236  Real gs = g.apply(s);
237  Real modelVal = fold - pRed;
238  Real theta = (one-eta2)*gs/((one-eta2)*(fold+gs)+eta2*modelVal-ftrial);
239  if ( print ) {
240  outStream << " In TrustRegionUtilities::interpolateRadius" << std::endl;
241  outStream << " Interpolation model value: " << modelVal << std::endl;
242  outStream << " Interpolation step length: " << theta << std::endl;
243  outStream << std::endl;
244  }
245  return std::min(gamma1*std::min(snorm,del),std::max(gamma0,theta)*del);
246 }
247 
248 } // namespace TrustRegion
249 } // namespace ROL
250 
251 
252 #endif
Provides the interface to evaluate objective functions.
Real initialRadius(int &nfval, const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &Bg, const Real fx, const Real gnorm, Objective< Real > &obj, TrustRegionModel_U< Real > &model, const Real delMax, std::ostream &outStream, const bool print=false)
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Real interpolateRadius(const Vector< Real > &g, const Vector< Real > &s, const Real snorm, const Real pRed, const Real fold, const Real ftrial, const Real del, const Real gamma0, const Real gamma1, const Real eta2, std::ostream &outStream=std::cout, const bool print=false)
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:238
virtual void setData(Objective< Real > &obj, const Vector< Real > &x, const Vector< Real > &g)
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
ETRFlag
Enumation of flags used by trust-region solvers.
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
Provides the interface to evaluate trust-region model functions.
std::string ETRFlagToString(ETRFlag trf)
void analyzeRatio(Real &rho, ETRFlag &flag, const Real fold, const Real ftrial, const Real pRed, const Real epsi, std::ostream &outStream=std::cout, const bool print=false)