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