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