ROL
ROL_BPOEObjective.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_BPOEOBJECTIVE_HPP
45 #define ROL_BPOEOBJECTIVE_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 BPOEObjective : public Objective<Real> {
57 private:
58  Teuchos::RCP<ParametrizedObjective<Real> > ParametrizedObjective_;
59 
60  Real order_;
61  Real threshold_;
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 
79  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 public:
148  virtual ~BPOEObjective() {}
149 
150  BPOEObjective( const Teuchos::RCP<ParametrizedObjective<Real> > &pObj,
151  const Real order, const Real threshold,
152  const Teuchos::RCP<SampleGenerator<Real> > &vsampler,
153  const Teuchos::RCP<SampleGenerator<Real> > &gsampler,
154  const Teuchos::RCP<SampleGenerator<Real> > &hsampler,
155  const bool storage = true )
156  : ParametrizedObjective_(pObj), order_(order), threshold_(threshold),
157  ValueSampler_(vsampler), GradientSampler_(gsampler), HessianSampler_(hsampler),
158  storage_(storage), initialized_(false) {
159  value_storage_.clear();
160  gradient_storage_.clear();
161  }
162 
163  BPOEObjective( const Teuchos::RCP<ParametrizedObjective<Real> > &pObj,
164  const Real order, const Real threshold,
165  const Teuchos::RCP<SampleGenerator<Real> > &vsampler,
166  const Teuchos::RCP<SampleGenerator<Real> > &gsampler,
167  const bool storage = true )
168  : ParametrizedObjective_(pObj), order_(order), threshold_(threshold),
169  ValueSampler_(vsampler), GradientSampler_(gsampler), HessianSampler_(gsampler),
170  storage_(storage), initialized_(false) {
171  value_storage_.clear();
172  gradient_storage_.clear();
173  }
174 
175  BPOEObjective( const Teuchos::RCP<ParametrizedObjective<Real> > &pObj,
176  const Real order, const Real threshold,
177  const Teuchos::RCP<SampleGenerator<Real> > &sampler,
178  const bool storage = true )
179  : ParametrizedObjective_(pObj), order_(order), threshold_(threshold),
180  ValueSampler_(sampler), GradientSampler_(sampler), HessianSampler_(sampler),
181  storage_(storage), initialized_(false) {
182  value_storage_.clear();
183  gradient_storage_.clear();
184  }
185 
186  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
187  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
188  unwrap_const_CVaR_vector(xvec,xvar,x);
189  ParametrizedObjective_->update(*xvec,flag,iter);
190  ValueSampler_->update(*xvec);
191  if ( storage_ ) {
192  value_storage_.clear();
193  }
194  if ( flag ) {
195  GradientSampler_->update(*xvec);
196  HessianSampler_->update(*xvec);
197  if ( storage_ ) {
198  gradient_storage_.clear();
199  }
200  }
201  }
202 
203  Real value( const Vector<Real> &x, Real &tol ) {
204  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
205  unwrap_const_CVaR_vector(xvec,xvar,x);
206  // Initialize storage
207  std::vector<Real> point;
208  Real weight = 0.0, myval = 0.0, pval = 0.0, val = 0.0, bp = 0.0;
209  int start = ValueSampler_->start(), end = ValueSampler_->numMySamples();
210  for ( int i = start; i < end; i++ ) {
211  weight = ValueSampler_->getMyWeight(i);
212  point = ValueSampler_->getMyPoint(i);
213  // Compute f(xvec,xi)
214  getValue(pval,*xvec,point,tol);
215  bp = xvar*(pval-threshold_)+1.0;
216  if ( bp > 0.0 ) {
217  // Build partial sum depending on value
218  myval += weight*((order_==1.0) ? bp
219  : std::pow(bp,order_));
220  }
221  }
222  // Update expected value
223  ValueSampler_->sumAll(&myval,&val,1);
224  // Return BPOE value
225  return ((order_==1.0) ? val : std::pow(val,1.0/order_));
226  }
227 
228  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
229  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
230  unwrap_const_CVaR_vector(xvec,xvar,x);
231  RiskVector<Real> &gc = Teuchos::dyn_cast<RiskVector<Real> >(
232  Teuchos::dyn_cast<Vector<Real> >(g));
233  // Initialize storage
234  g.zero(); sumGrad0_->zero(); pointGrad_->zero();
235  std::vector<Real> point, val(2,0.0), myval(2,0.0);
236  Real weight = 0.0, pval = 0.0, pvalp0 = 0.0, pvalp1 = 0.0, bp = 0.0;
237  int start = GradientSampler_->start(), end = GradientSampler_->numMySamples();
238  for ( int i = start; i < end; i++ ) {
239  weight = GradientSampler_->getMyWeight(i);
240  point = GradientSampler_->getMyPoint(i);
241  // Compute the value of f(xvec,xi)
242  getValue(pval,*xvec,point,tol);
243  bp = xvar*(pval-threshold_)+1.0;
244  if ( bp > 0.0 ) {
245  // Compute max(0,f(xvec,xi)-xvar)^order
246  pvalp0 = ((order_==1.0) ? bp
247  : std::pow(bp,order_));
248  pvalp1 = ((order_==1.0) ? 1.0
249  : ((order_==2.0) ? bp
250  : std::pow(bp,order_-1.0)));
251  // Build partial sums depending on value
252  myval[0] += weight*pvalp0;
253  myval[1] += weight*pvalp1*(pval-threshold_);
254  // Compute gradient of f(xvec,xi)
255  getGradient(*pointGrad_,*xvec,point,tol);
256  // Build partial sum depending on gradient
257  sumGrad0_->axpy(weight*pvalp1,*pointGrad_);
258  }
259  }
260  // Combine partial sums
261  GradientSampler_->sumAll(&myval[0],&val[0],2);
262  // Compute VaR gradient and BPOE gradient
263  Real gvar = 0.0; gradient0_->zero();
264  if ( std::abs(val[0]) >= ROL_EPSILON) {
265  GradientSampler_->sumAll(*sumGrad0_,*gradient0_);
266  Real norm = std::pow(val[0],(order_-1.0)/order_);
267  gradient0_->scale(xvar/norm);
268  gvar = val[1]/norm;
269  }
270  // Set gradient components of CVaR vector
271  gc.setStatistic(gvar);
272  gc.setVector(*(Teuchos::rcp_dynamic_cast<Vector<Real> >(gradient0_)));
273  }
274 
275  void hessVec( Vector<Real> &hv, const Vector<Real> &v,
276  const Vector<Real> &x, Real &tol ) {
277  Teuchos::RCP<Vector<Real> > xvec; Real xvar = 0.0;
278  unwrap_const_CVaR_vector(xvec,xvar,x);
279  Teuchos::RCP<Vector<Real> > vvec; Real vvar = 0.0;
280  unwrap_const_CVaR_vector(vvec,vvar,v);
281  RiskVector<Real> &hvc = Teuchos::dyn_cast<RiskVector<Real> >(
282  Teuchos::dyn_cast<Vector<Real> >(hv));
283  // Initialize storage
284  hv.zero(); sumHess_->zero(); hessvec_->zero();
285  sumGrad0_->zero(); sumGrad1_->zero(); sumGrad2_->zero();
286  gradient0_->zero(); gradient1_->zero(); gradient2_->zero();
287  std::vector<Real> point, val(5,0.0), myval(5,0.0);
288  Real pval = 0.0, pvalp0 = 0.0, pvalp1 = 0.0, pvalp2 = 0.0;
289  Real weight = 0.0, gv = 0.0, bp = 0.0;
290  int start = HessianSampler_->start(), end = HessianSampler_->numMySamples();
291  for ( int i = start; i < end; i++ ) {
292  // Get sample and associated probability
293  weight = HessianSampler_->getMyWeight(i);
294  point = HessianSampler_->getMyPoint(i);
295  // Compute the value of f(xvec,xi)
296  getValue(pval,*xvec,point,tol);
297  bp = xvar*(pval-threshold_)+1.0;
298  if ( bp > 0.0 ) {
299  // Compute max(0,f(xvec,xi)-xvar)^order
300  pvalp0 = ((order_==1.0) ? bp
301  : std::pow(bp,order_));
302  pvalp1 = ((order_==1.0) ? 1.0
303  : ((order_==2.0) ? bp
304  : std::pow(bp,order_-1.0)));
305  pvalp2 = ((order_==1.0) ? 0.0
306  : ((order_==2.0) ? 1.0
307  : ((order_==3.0) ? bp
308  : std::pow(bp,order_-2.0))));
309  // Build partial sums depending on value
310  myval[0] += weight*pvalp0;
311  myval[1] += weight*pvalp1*(pval-threshold_);
312  myval[2] += weight*pvalp2*(pval-threshold_)*(pval-threshold_);
313  // Compute the gradient and directional derivative of f(xvec,xi)
314  getGradient(*pointGrad_,*xvec,point,tol);
315  gv = pointGrad_->dot(vvec->dual());
316  // Build partial sums depending on gradient
317  myval[3] += weight*pvalp1*gv;
318  myval[4] += weight*pvalp2*(pval-threshold_)*gv;
319  sumGrad0_->axpy(weight*pvalp1,*pointGrad_);
320  sumGrad1_->axpy(weight*pvalp2*(pval-threshold_),*pointGrad_);
321  sumGrad2_->axpy(weight*pvalp2*gv,*pointGrad_);
322  // Compute the hessian of f(xvec,xi) in the direction vvec
323  getHessVec(*pointHess_,*vvec,*xvec,point,tol);
324  // Build partial sum depending on the hessian
325  sumHess_->axpy(weight*pvalp1,*pointHess_);
326  }
327  }
328  // Compile partial sums
329  HessianSampler_->sumAll(&myval[0],&val[0],5);
330  Real hvar = 0.0; hessvec_->zero();
331  if ( std::abs(val[0]) >= ROL_EPSILON ) {
332  HessianSampler_->sumAll(*sumGrad0_,*gradient0_);
333  HessianSampler_->sumAll(*sumGrad1_,*gradient1_);
334  HessianSampler_->sumAll(*sumGrad2_,*gradient2_);
335  HessianSampler_->sumAll(*sumHess_,*hessvec_);
336  // Compute VaR Hessian-times-a-vector and BPOE Hessian-times-a-vector
337  Real norm0 = ((order_==1.0) ? 1.0
338  : ((order_==2.0) ? std::sqrt(val[0])
339  : std::pow(val[0],(order_-1.0)/order_)));
340  Real norm1 = ((order_==1.0) ? val[0]
341  : std::pow(val[0],(2.0*order_-1.0)/order_));
342  hvar = (order_-1.0)*((val[2]/norm0 - val[1]*val[1]/norm1)*vvar
343  +xvar*(val[4]/norm0 - val[3]*val[1]/norm1))
344  +(val[3]/norm0);
345  hessvec_->scale(xvar/norm0); //(order_-1.0)/norm0);
346  Real coeff = -(order_-1.0)*xvar*(xvar*val[3]+vvar*val[1])/norm1+vvar/norm0;
347  hessvec_->axpy(coeff,*gradient0_);
348  hessvec_->axpy((order_-1.0)*vvar*xvar/norm0,*gradient1_);
349  hessvec_->axpy((order_-1.0)*xvar*xvar/norm0,*gradient2_);
350  }
351  // Set gradient components of CVaR vector
352  hvc.setStatistic(hvar);
353  hvc.setVector(*(Teuchos::rcp_dynamic_cast<Vector<Real> >(hessvec_)));
354  }
355 
356  virtual void precond( Vector<Real> &Pv, const Vector<Real> &v,
357  const Vector<Real> &x, Real &tol ) {
358  Pv.set(v.dual());
359  }
360 };
361 
362 }
363 
364 #endif
Provides the interface to evaluate objective functions.
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 > > gradient1_
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)
Teuchos::RCP< SampleGenerator< Real > > HessianSampler_
Teuchos::RCP< Vector< Real > > pointGrad_
Teuchos::RCP< Vector< Real > > sumGrad1_
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
void initialize(const Vector< Real > &x)
void getGradient(Vector< Real > &g, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
Teuchos::RCP< Vector< Real > > sumGrad0_
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.
void getValue(Real &val, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
Teuchos::RCP< Vector< Real > > pointHess_
std::map< std::vector< Real >, Real > value_storage_
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
std::map< std::vector< Real >, Teuchos::RCP< Vector< Real > > > gradient_storage_
void setVector(const Vector< Real > &vec)
Teuchos::RCP< const Vector< Real > > getVector() const
Teuchos::RCP< SampleGenerator< Real > > ValueSampler_
BPOEObjective(const Teuchos::RCP< ParametrizedObjective< Real > > &pObj, const Real order, const Real threshold, const Teuchos::RCP< SampleGenerator< Real > > &vsampler, const Teuchos::RCP< SampleGenerator< Real > > &gsampler, const Teuchos::RCP< SampleGenerator< Real > > &hsampler, const bool storage=true)
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.
BPOEObjective(const Teuchos::RCP< ParametrizedObjective< Real > > &pObj, const Real order, const Real threshold, const Teuchos::RCP< SampleGenerator< Real > > &vsampler, const Teuchos::RCP< SampleGenerator< Real > > &gsampler, const bool storage=true)
Teuchos::RCP< Vector< Real > > gradient2_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Teuchos::RCP< Vector< Real > > sumHess_
Teuchos::RCP< Vector< Real > > sumGrad2_
Teuchos::RCP< Vector< Real > > gradient0_
BPOEObjective(const Teuchos::RCP< ParametrizedObjective< Real > > &pObj, const Real order, const Real threshold, const Teuchos::RCP< SampleGenerator< Real > > &sampler, const bool storage=true)
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void getHessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, const std::vector< Real > &param, Real &tol)
Teuchos::RCP< SampleGenerator< Real > > GradientSampler_
Teuchos::RCP< Vector< Real > > hessvec_
Teuchos::RCP< ParametrizedObjective< Real > > ParametrizedObjective_
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118