ROL
ROL_HMCRObjective.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_HMCROBJECTIVE_HPP
45 #define ROL_HMCROBJECTIVE_HPP
46 
47 #include "Teuchos_RCP.hpp"
48 #include "ROL_RiskVector.hpp"
49 #include "ROL_Objective.hpp"
51 #include "ROL_SampleGenerator.hpp"
52 
53 namespace ROL {
54 
55 template<class Real>
56 class HMCRObjective : public Objective<Real> {
57 private:
58  Teuchos::RCP<ParametrizedObjective<Real> > ParametrizedObjective_;
59 
60  Real order_;
61  Real prob_;
62 
63  Teuchos::RCP<SampleGenerator<Real> > ValueSampler_;
64  Teuchos::RCP<SampleGenerator<Real> > GradientSampler_;
65  Teuchos::RCP<SampleGenerator<Real> > HessianSampler_;
66 
67  Teuchos::RCP<Vector<Real> > pointGrad_;
68  Teuchos::RCP<Vector<Real> > pointHess_;
69 
70  Teuchos::RCP<Vector<Real> > gradient0_;
71  Teuchos::RCP<Vector<Real> > sumGrad0_;
72  Teuchos::RCP<Vector<Real> > gradient1_;
73  Teuchos::RCP<Vector<Real> > sumGrad1_;
74  Teuchos::RCP<Vector<Real> > gradient2_;
75  Teuchos::RCP<Vector<Real> > sumGrad2_;
76  Teuchos::RCP<Vector<Real> > hessvec_;
77  Teuchos::RCP<Vector<Real> > sumHess_;
78 
80  bool storage_;
81 
82  std::map<std::vector<Real>,Real> value_storage_;
83  std::map<std::vector<Real>,Teuchos::RCP<Vector<Real> > > gradient_storage_;
84 
85  void initialize(const Vector<Real> &x) {
86  pointGrad_ = x.dual().clone();
87  pointHess_ = x.dual().clone();
88  gradient0_ = x.dual().clone();
89  sumGrad0_ = x.dual().clone();
90  gradient1_ = x.dual().clone();
91  sumGrad1_ = x.dual().clone();
92  gradient2_ = x.dual().clone();
93  sumGrad2_ = x.dual().clone();
94  hessvec_ = x.dual().clone();
95  sumHess_ = x.dual().clone();
96  initialized_ = true;
97  }
98 
99  void unwrap_const_CVaR_vector(Teuchos::RCP<Vector<Real> > &xvec, Real &xvar,
100  const Vector<Real> &x) {
101  xvec = Teuchos::rcp_const_cast<Vector<Real> >(
102  Teuchos::dyn_cast<const RiskVector<Real> >(
103  Teuchos::dyn_cast<const Vector<Real> >(x)).getVector());
104  xvar = Teuchos::dyn_cast<const RiskVector<Real> >(
105  Teuchos::dyn_cast<const Vector<Real> >(x)).getStatistic();
106  if ( !initialized_ ) {
107  initialize(*xvec);
108  }
109  }
110 
111  void getValue(Real &val, const Vector<Real> &x,
112  const std::vector<Real> &param, Real &tol) {
113  if ( storage_ && value_storage_.count(param) ) {
114  val = value_storage_[param];
115  }
116  else {
117  ParametrizedObjective_->setParameter(param);
118  val = ParametrizedObjective_->value(x,tol);
119  if ( storage_ ) {
120  value_storage_.insert(std::pair<std::vector<Real>,Real>(param,val));
121  }
122  }
123  }
124 
126  const std::vector<Real> &param, Real &tol) {
127  if ( storage_ && gradient_storage_.count(param) ) {
128  g.set(*(gradient_storage_[param]));
129  }
130  else {
131  ParametrizedObjective_->setParameter(param);
132  ParametrizedObjective_->gradient(g,x,tol);
133  if ( storage_ ) {
134  Teuchos::RCP<Vector<Real> > tmp = g.clone();
135  gradient_storage_.insert(std::pair<std::vector<Real>,Teuchos::RCP<Vector<Real> > >(param,tmp));
136  gradient_storage_[param]->set(g);
137  }
138  }
139  }
140 
141  void getHessVec(Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x,
142  const std::vector<Real> &param, Real &tol) {
143  ParametrizedObjective_->setParameter(param);
144  ParametrizedObjective_->hessVec(hv,v,x,tol);
145  }
146 
147 
148 public:
149  virtual ~HMCRObjective() {}
150 
152  Real order, Real prob,
153  Teuchos::RCP<SampleGenerator<Real> > &vsampler,
154  Teuchos::RCP<SampleGenerator<Real> > &gsampler,
155  Teuchos::RCP<SampleGenerator<Real> > &hsampler,
156  bool storage = true )
157  : ParametrizedObjective_(pObj),
158  ValueSampler_(vsampler), GradientSampler_(gsampler), HessianSampler_(hsampler),
159  initialized_(false), storage_(storage) {
160  order_ = ((order < 1.0) ? 1.0 : order);
161  prob_ = ((prob > 1.0) ? 1.0 : ((prob < 0.0) ? 0.0 : prob));
162  value_storage_.clear();
163  gradient_storage_.clear();
164  }
165 
167  Real order, Real prob,
168  Teuchos::RCP<SampleGenerator<Real> > &vsampler,
169  Teuchos::RCP<SampleGenerator<Real> > &gsampler,
170  bool storage = true )
171  : ParametrizedObjective_(pObj),
172  ValueSampler_(vsampler), GradientSampler_(gsampler), HessianSampler_(gsampler),
173  initialized_(false), storage_(storage) {
174  order_ = ((order < 1.0) ? 1.0 : order);
175  prob_ = ((prob > 1.0) ? 1.0 : ((prob < 0.0) ? 0.0 : prob));
176  value_storage_.clear();
177  gradient_storage_.clear();
178  }
179 
181  Real order, Real prob,
182  Teuchos::RCP<SampleGenerator<Real> > &sampler,
183  bool storage = true )
184  : ParametrizedObjective_(pObj), order_(order), prob_(prob),
185  ValueSampler_(sampler), GradientSampler_(sampler), HessianSampler_(sampler),
186  initialized_(false), storage_(storage) {
187  order_ = ((order < 1.0) ? 1.0 : order);
188  prob_ = ((prob > 1.0) ? 1.0 : ((prob < 0.0) ? 0.0 : prob));
189  value_storage_.clear();
190  gradient_storage_.clear();
191  }
192 
193  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
194  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
195  unwrap_const_CVaR_vector(xvec,xvar,x);
196  ParametrizedObjective_->update(*xvec,flag,iter);
197  ValueSampler_->update(*xvec);
198  if ( storage_ ) {
199  value_storage_.clear();
200  }
201  if ( flag ) {
202  GradientSampler_->update(*xvec);
203  HessianSampler_->update(*xvec);
204  if ( storage_ ) {
205  gradient_storage_.clear();
206  }
207  }
208  }
209 
210  Real value( const Vector<Real> &x, Real &tol ) {
211  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
212  unwrap_const_CVaR_vector(xvec,xvar,x);
213  // Initialize storage
214  std::vector<Real> point;
215  Real weight = 0.0, myval = 0.0, pval = 0.0, val = 0.0;
216  int start = ValueSampler_->start(), end = ValueSampler_->numMySamples();
217  for ( int i = start; i < end; i++ ) {
218  weight = ValueSampler_->getMyWeight(i);
219  point = ValueSampler_->getMyPoint(i);
220  // Compute f(xvec,xi)
221  getValue(pval,*xvec,point,tol);
222  if ( pval > xvar ) {
223  // Build partial sum depending on value
224  myval += weight*((order_ == 1.0) ? pval-xvar
225  : std::pow(pval-xvar,order_));
226  }
227  }
228  // Update expected value
229  ValueSampler_->sumAll(&myval,&val,1);
230  // Return HMCR value
231  if (std::abs(val) < ROL_EPSILON) {
232  return xvar;
233  }
234  return xvar + ((order_ == 1.0) ? val
235  : ((order_ == 2.0) ? std::sqrt(val)
236  : std::pow(val,1.0/order_)))/(1.0-prob_);
237  }
238 
239  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
240  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
241  unwrap_const_CVaR_vector(xvec,xvar,x);
242  RiskVector<Real> &gc = Teuchos::dyn_cast<RiskVector<Real> >(
243  Teuchos::dyn_cast<Vector<Real> >(g));
244  // Initialize storage
245  g.zero(); sumGrad0_->zero();
246  std::vector<Real> point, val(2,0.0), myval(2,0.0);
247  Real weight = 0.0, pval = 0.0, pvalp0 = 0.0, pvalp1 = 0.0;
248  int start = GradientSampler_->start(), end = GradientSampler_->numMySamples();
249  for ( int i = start; i < end; i++ ) {
250  weight = GradientSampler_->getMyWeight(i);
251  point = GradientSampler_->getMyPoint(i);
252  // Compute the value of f(xvec,xi)
253  getValue(pval,*xvec,point,tol);
254  if ( pval > xvar ) {
255  // Compute max(0,f(xvec,xi)-xvar)^order
256  pvalp0 = ((order_ == 1.0) ? pval-xvar
257  : std::pow(pval-xvar,order_));
258  pvalp1 = ((order_ == 1.0) ? 1.0
259  : ((order_ == 2.0) ? pval-xvar
260  : std::pow(pval-xvar,order_-1.0)));
261  // Build partial sums depending on value
262  myval[0] += weight*pvalp0;
263  myval[1] += weight*pvalp1;
264  // Compute gradient of f(xvec,xi)
265  getGradient(*pointGrad_,*xvec,point,tol);
266  // Build partial sum depending on gradient
267  sumGrad0_->axpy(weight*pvalp1,*pointGrad_);
268  }
269  }
270  Real gvar = 1.0; gradient0_->zero();
271  // Combine partial sums
272  GradientSampler_->sumAll(&myval[0],&val[0],2);
273  if (std::abs(val[0]) >= ROL_EPSILON) {
274  GradientSampler_->sumAll(*sumGrad0_,*gradient0_);
275  // Compute VaR gradient and HMCR gradient
276  Real norm = ((order_ == 1.0) ? 1.0
277  : ((order_ == 2.0) ? std::sqrt(val[0])
278  : std::pow(val[0],(order_-1.0)/order_)));
279  gvar -= val[1]/((1.0-prob_)*norm);
280  gradient0_->scale(1.0/((1.0-prob_)*norm));
281  }
282  // Set gradient components of CVaR vector
283  gc.setStatistic(gvar);
284  gc.setVector(*(Teuchos::rcp_dynamic_cast<Vector<Real> >(gradient0_)));
285  }
286 
287  void hessVec( Vector<Real> &hv, const Vector<Real> &v,
288  const Vector<Real> &x, Real &tol ) {
289  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
290  unwrap_const_CVaR_vector(xvec,xvar,x);
291  Teuchos::RCP<Vector<Real> > vvec; Real vvar = 0.0;
292  unwrap_const_CVaR_vector(vvec,vvar,v);
293  RiskVector<Real> &hvc = Teuchos::dyn_cast<RiskVector<Real> >(
294  Teuchos::dyn_cast<Vector<Real> >(hv));
295  // Initialize storage
296  hv.zero();
297  sumGrad0_->zero(); sumGrad1_->zero(); sumGrad2_->zero(); sumHess_->zero();
298  gradient0_->zero(); gradient1_->zero(); gradient2_->zero();
299  Real weight = 0.0;
300  std::vector<Real> point;
301  Real pval = 0.0, pvalp0 = 0.0, pvalp1 = 0.0, pvalp2 = 0.0, gv = 0.0;
302  std::vector<Real> val(5,0.0), myval(5,0.0);
303  int start = HessianSampler_->start(), end = HessianSampler_->numMySamples();
304  for ( int i = start; i < end; i++ ) {
305  // Get sample and associated probability
306  weight = HessianSampler_->getMyWeight(i);
307  point = HessianSampler_->getMyPoint(i);
308  // Compute the value of f(xvec,xi)
309  getValue(pval,*xvec,point,tol);
310  if ( pval > xvar ) {
311  // Compute max(0,f(xvec,xi)-xvar)^order
312  pvalp0 = ((order_ == 1.0) ? pval-xvar
313  : std::pow(pval-xvar,order_));
314  pvalp1 = ((order_ == 1.0) ? 1.0
315  : ((order_ == 2.0) ? pval-xvar
316  : std::pow(pval-xvar,order_-1.0)));
317  pvalp2 = ((order_ == 1.0) ? 0.0
318  : ((order_ == 2.0) ? 1.0
319  : ((order_ == 3.0) ? pval-xvar
320  : std::pow(pval-xvar,order_-2.0))));
321  // Build partial sums depending on value
322  myval[0] += weight*pvalp0;
323  myval[1] += weight*pvalp1;
324  myval[2] += weight*pvalp2;
325  // Compute the gradient and directional derivative of f(xvec,xi)
326  getGradient(*pointGrad_,*xvec,point,tol);
327  gv = pointGrad_->dot(vvec->dual());
328  // Build partial sums depending on gradient
329  myval[3] += weight*pvalp1*gv;
330  myval[4] += weight*pvalp2*gv;
331  sumGrad0_->axpy(weight*pvalp1,*pointGrad_);
332  sumGrad1_->axpy(weight*pvalp2,*pointGrad_);
333  sumGrad2_->axpy(weight*pvalp2*gv,*pointGrad_);
334  // Compute the hessian of f(xvec,xi) in the direction vvec
335  getHessVec(*pointHess_,*vvec,*xvec,point,tol);
336  // Build partial sum depending on the hessian
337  sumHess_->axpy(weight*pvalp1,*pointHess_);
338  }
339  }
340  Real hvar = 0.0; hessvec_->zero();
341  HessianSampler_->sumAll(&myval[0],&val[0],5);
342  if (std::abs(val[0]) >= ROL_EPSILON) {
343  // Compile partial sums
344  HessianSampler_->sumAll(*sumGrad0_,*gradient0_);
345  HessianSampler_->sumAll(*sumGrad1_,*gradient1_);
346  HessianSampler_->sumAll(*sumGrad2_,*gradient2_);
347  HessianSampler_->sumAll(*sumHess_,*hessvec_);
348  // Compute VaR Hessian-times-a-vector and HMCR Hessian-times-a-vector
349  Real norm0 = (1.0-prob_)*((order_ == 1.0) ? 1.0
350  : ((order_ == 2.0) ? std::sqrt(val[0])
351  : std::pow(val[0],(order_-1.0)/order_)));
352  Real norm1 = (1.0-prob_)*((order_ == 1.0) ? 1.0
353  : std::pow(val[0],(2.0*order_-1.0)/order_));
354  hvar = (order_-1.0)*((val[2]/norm0 - val[1]*val[1]/norm1)*vvar
355  -(val[4]/norm0 - val[3]*val[1]/norm1));
356  hessvec_->scale(1.0/norm0); //(order_-1.0)/norm0);
357  hessvec_->axpy(-(order_-1.0)*vvar/norm0,*gradient1_);
358  hessvec_->axpy((order_-1.0)*(vvar*val[1]-val[3])/norm1,*gradient0_);
359  hessvec_->axpy((order_-1.0)/norm0,*gradient2_);
360  }
361  // Set gradient components of CVaR vector
362  hvc.setStatistic(hvar);
363  hvc.setVector(*(Teuchos::rcp_dynamic_cast<Vector<Real> >(hessvec_)));
364  }
365 
366  virtual void precond( Vector<Real> &Pv, const Vector<Real> &v,
367  const Vector<Real> &x, Real &tol ) {
368  Pv.set(v.dual());
369  }
370 };
371 
372 }
373 
374 #endif
void getValue(Real &val, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
Provides the interface to evaluate objective functions.
Teuchos::RCP< SampleGenerator< Real > > GradientSampler_
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< ParametrizedObjective< Real > > ParametrizedObjective_
Teuchos::RCP< SampleGenerator< Real > > ValueSampler_
void getHessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
Teuchos::RCP< Vector< Real > > sumGrad1_
Teuchos::RCP< Vector< Real > > gradient1_
Teuchos::RCP< Vector< Real > > hessvec_
Teuchos::RCP< Vector< Real > > sumGrad2_
Teuchos::RCP< Vector< Real > > gradient0_
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
const Real getStatistic(const int i=0) const
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< Vector< Real > > pointGrad_
std::map< std::vector< Real >, Real > value_storage_
Teuchos::RCP< SampleGenerator< Real > > HessianSampler_
std::map< std::vector< Real >, Teuchos::RCP< Vector< Real > > > gradient_storage_
void setVector(const Vector< Real > &vec)
Teuchos::RCP< const Vector< Real > > getVector() const
void initialize(const Vector< Real > &x)
void setStatistic(const Real stat)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< Vector< Real > > sumHess_
HMCRObjective(Teuchos::RCP< ParametrizedObjective< Real > > &pObj, Real order, Real prob, Teuchos::RCP< SampleGenerator< Real > > &sampler, bool storage=true)
Teuchos::RCP< Vector< Real > > pointHess_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
HMCRObjective(Teuchos::RCP< ParametrizedObjective< Real > > &pObj, Real order, Real prob, Teuchos::RCP< SampleGenerator< Real > > &vsampler, Teuchos::RCP< SampleGenerator< Real > > &gsampler, bool storage=true)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Teuchos::RCP< Vector< Real > > gradient2_
void getGradient(Vector< Real > &g, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void unwrap_const_CVaR_vector(Teuchos::RCP< Vector< Real > > &xvec, Real &xvar, const Vector< Real > &x)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Teuchos::RCP< Vector< Real > > sumGrad0_
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118
HMCRObjective(Teuchos::RCP< ParametrizedObjective< Real > > &pObj, Real order, Real prob, Teuchos::RCP< SampleGenerator< Real > > &vsampler, Teuchos::RCP< SampleGenerator< Real > > &gsampler, Teuchos::RCP< SampleGenerator< Real > > &hsampler, bool storage=true)