ROL
ROL_HMCR.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_HMCR_HPP
45#define ROL_HMCR_HPP
46
48#include "ROL_PlusFunction.hpp"
49#include "ROL_RiskVector.hpp"
50
72namespace ROL {
73
74template<class Real>
75class HMCR : public RandVarFunctional<Real> {
76private:
77 // User inputs
78 ROL::Ptr<PlusFunction<Real> > plusFunction_;
79 Real prob_;
80 Real lambda_;
81 unsigned order_;
82
83 // 1/(1-prob)
84 Real coeff_;
85
86 // Temporary vector storage
87 ROL::Ptr<Vector<Real> > mDualVector_;
88
89 // Temporary scalar storage
90 Real pnorm_;
91 Real coeff0_;
92 Real coeff1_;
93 Real coeff2_;
94
95 // Flag to initialized vector storage
97
98 using RandVarFunctional<Real>::val_;
99 using RandVarFunctional<Real>::gv_;
100 using RandVarFunctional<Real>::g_;
101 using RandVarFunctional<Real>::hv_;
103
104 using RandVarFunctional<Real>::point_;
105 using RandVarFunctional<Real>::weight_;
106
111
112 void checkInputs(void) const {
113 const Real zero(0), one(1);
114 ROL_TEST_FOR_EXCEPTION((prob_ <= zero) || (prob_ >= one), std::invalid_argument,
115 ">>> ERROR (ROL::HMCR): Confidence level must be between 0 and 1!");
116 ROL_TEST_FOR_EXCEPTION((lambda_ < zero) || (lambda_ > one), std::invalid_argument,
117 ">>> ERROR (ROL::HMCR): Convex combination parameter must be positive!");
118 ROL_TEST_FOR_EXCEPTION((order_ < 2), std::invalid_argument,
119 ">>> ERROR (ROL::HMCR): Norm order is less than 2!");
120 ROL_TEST_FOR_EXCEPTION(plusFunction_ == ROL::nullPtr, std::invalid_argument,
121 ">>> ERROR (ROL::HMCR): PlusFunction pointer is null!");
122 }
123
124public:
134 HMCR( const Real prob, const Real lambda, const unsigned order,
135 const ROL::Ptr<PlusFunction<Real> > &pf )
136 : RandVarFunctional<Real>(),
137 plusFunction_(pf), prob_(prob), lambda_(lambda), order_(order),
138 pnorm_(0), coeff0_(0), coeff1_(0), coeff2_(0),
139 HMCR_firstReset_(true) {
140 checkInputs();
141 const Real one(1);
142 coeff_ = one/(one-prob_);
143 }
144
156 HMCR( ROL::ParameterList &parlist )
157 : RandVarFunctional<Real>(),
158 pnorm_(0), coeff0_(0), coeff1_(0), coeff2_(0),
159 HMCR_firstReset_(true) {
160 ROL::ParameterList &list
161 = parlist.sublist("SOL").sublist("Risk Measure").sublist("HMCR");
162 // Check HMCR inputs
163 prob_ = list.get<Real>("Confidence Level");
164 lambda_ = list.get<Real>("Convex Combination Parameter");
165 order_ = (unsigned)list.get<int>("Order",2);
166 // Build (approximate) plus function
167 plusFunction_ = ROL::makePtr<PlusFunction<Real> >(list);
168 // Check inputs
169 checkInputs();
170 const Real one(1);
171 coeff_ = one/(one-prob_);
172 }
173
174 void initialize(const Vector<Real> &x) {
176 // Initialize additional vector storage
177 if ( HMCR_firstReset_ ) {
178 mDualVector_ = x.dual().clone();
179 HMCR_firstReset_ = false;
180 }
181 // Zero temporary storage
182 const Real zero(0);
183 mDualVector_->zero();
184 pnorm_ = zero; coeff0_ = zero;
186 }
187
189 const Vector<Real> &x,
190 const std::vector<Real> &xstat,
191 Real &tol) {
192 const Real rorder = static_cast<Real>(order_);
193 // Expected value
194 Real val = computeValue(obj,x,tol);
195 val_ += weight_*val;
196 // Higher moment
197 Real pf = plusFunction_->evaluate(val-xstat[0],0);
198 pnorm_ += weight_*std::pow(pf,rorder);
199 }
200
201 Real getValue(const Vector<Real> &x,
202 const std::vector<Real> &xstat,
203 SampleGenerator<Real> &sampler) {
204 const Real one(1);
205 const Real power = one/static_cast<Real>(order_);
206 std::vector<Real> val_in(2), val_out(2);
207 val_in[0] = val_;
208 val_in[1] = pnorm_;
209 sampler.sumAll(&val_in[0],&val_out[0],2);
210 return (one-lambda_)*val_out[0]
211 + lambda_*(xstat[0] + coeff_*std::pow(val_out[1],power));
212 }
213
215 const Vector<Real> &x,
216 const std::vector<Real> &xstat,
217 Real &tol) {
218 const Real one(1);
219 const Real rorder0 = static_cast<Real>(order_);
220 const Real rorder1 = rorder0 - one;
221 // Expected value
222 computeGradient(*dualVector_,obj,x,tol);
223 g_->axpy(weight_,*dualVector_);
224 // Higher moment
225 Real val = computeValue(obj,x,tol);
226 Real pf0 = plusFunction_->evaluate(val-xstat[0],0);
227 Real pf1 = plusFunction_->evaluate(val-xstat[0],1);
228
229 Real pf0p0 = std::pow(pf0,rorder0);
230 Real pf0p1 = std::pow(pf0,rorder1);
231
232 pnorm_ += weight_*pf0p0;
233 coeff0_ += weight_*pf0p1*pf1;
234
235 hv_->axpy(weight_*pf0p1*pf1,*dualVector_);
236 }
237
239 std::vector<Real> &gstat,
240 const Vector<Real> &x,
241 const std::vector<Real> &xstat,
242 SampleGenerator<Real> &sampler) {
243 const Real zero(0), one(1);
244 std::vector<Real> val_in(2), val_out(2);
245 val_in[0] = pnorm_; val_in[1] = coeff0_;
246
247 sampler.sumAll(&val_in[0],&val_out[0],2);
248 sampler.sumAll(*g_,g);
249 g.scale(one-lambda_);
250 Real var = lambda_;
251 // If the higher moment term is positive then compute gradient
252 if ( val_out[0] > zero ) {
253 const Real rorder0 = static_cast<Real>(order_);
254 const Real rorder1 = rorder0 - one;
255 Real denom = std::pow(val_out[0],rorder1/rorder0);
256 // Sum higher moment contribution
257 dualVector_->zero();
258 sampler.sumAll(*hv_,*dualVector_);
259 g.axpy(lambda_*coeff_/denom,*dualVector_);
260 // Compute statistic gradient
261 var -= lambda_*coeff_*((denom > zero) ? val_out[1]/denom : zero);
262 }
263 gstat[0] = var;
264 }
265
267 const Vector<Real> &v,
268 const std::vector<Real> &vstat,
269 const Vector<Real> &x,
270 const std::vector<Real> &xstat,
271 Real &tol) {
272 const Real one(1);
273 const Real rorder0 = static_cast<Real>(order_);
274 const Real rorder1 = rorder0-one;
275 const Real rorder2 = rorder1-one;
276 // Expected value
277 computeHessVec(*dualVector_,obj,v,x,tol);
278 hv_->axpy(weight_,*dualVector_);
279 // Higher moment
280 Real val = computeValue(obj,x,tol);
281 Real pf0 = plusFunction_->evaluate(val-xstat[0],0);
282 Real pf1 = plusFunction_->evaluate(val-xstat[0],1);
283 Real pf2 = plusFunction_->evaluate(val-xstat[0],2);
284
285 Real pf0p0 = std::pow(pf0,rorder0);
286 Real pf0p1 = std::pow(pf0,rorder1);
287 Real pf0p2 = std::pow(pf0,rorder2);
288
289 Real scale1 = pf0p1*pf1;
290 g_->axpy(weight_*scale1,*dualVector_);
291
292 Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
293 Real scale0 = (rorder1*pf0p2*pf1*pf1 + pf0p1*pf2)*(gv-vstat[0]);
294
295 pnorm_ += weight_*pf0p0;
296 coeff0_ += weight_*scale0;
297 coeff1_ += weight_*scale1;
298 coeff2_ += weight_*rorder1*scale1*(vstat[0]-gv);
299
300 g_->axpy(weight_*scale0,*dualVector_);
301 mDualVector_->axpy(weight_*scale1,*dualVector_);
302 }
303
305 std::vector<Real> &hvstat,
306 const Vector<Real> &v,
307 const std::vector<Real> &vstat,
308 const Vector<Real> &x,
309 const std::vector<Real> &xstat,
310 SampleGenerator<Real> &sampler) {
311 const Real zero(0), one(1);
312 std::vector<Real> val_in(4), val_out(4);
313 val_in[0] = pnorm_; val_in[1] = coeff0_;
314 val_in[2] = coeff1_; val_in[3] = coeff2_;
315
316 sampler.sumAll(&val_in[0],&val_out[0],4);
317 sampler.sumAll(*hv_,hv);
318
319 Real var = zero;
320 hv.scale(one-lambda_);
321
322 if ( val_out[0] > zero ) {
323 const Real rorder0 = static_cast<Real>(order_);
324 const Real rorder1 = rorder0-one;
325 const Real rorder2 = rorder0 + rorder1;
326 const Real coeff = lambda_*coeff_;
327
328 Real denom1 = std::pow(val_out[0],rorder1/rorder0);
329 Real denom2 = std::pow(val_out[0],rorder2/rorder0);
330
331 dualVector_->zero();
332 sampler.sumAll(*g_,*dualVector_);
333 hv.axpy(coeff/denom1,*dualVector_);
334
335 dualVector_->zero();
337 hv.axpy(coeff*val_out[3]/denom2,*dualVector_);
338
339 var = -coeff*(val_out[1]/denom1 + val_out[3]*val_out[2]/denom2);
340 }
341 hvstat[0] = var;
342 }
343};
344
345}
346
347#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides an interface for a convex combination of the expected value and the higher moment coherent r...
Definition: ROL_HMCR.hpp:75
bool HMCR_firstReset_
Definition: ROL_HMCR.hpp:96
void initialize(const Vector< Real > &x)
Initialize temporary variables.
Definition: ROL_HMCR.hpp:174
HMCR(const Real prob, const Real lambda, const unsigned order, const ROL::Ptr< PlusFunction< Real > > &pf)
Constructor.
Definition: ROL_HMCR.hpp:134
ROL::Ptr< Vector< Real > > mDualVector_
Definition: ROL_HMCR.hpp:87
HMCR(ROL::ParameterList &parlist)
Constructor.
Definition: ROL_HMCR.hpp:156
Real coeff1_
Definition: ROL_HMCR.hpp:92
unsigned order_
Definition: ROL_HMCR.hpp:81
Real coeff_
Definition: ROL_HMCR.hpp:84
void getGradient(Vector< Real > &g, std::vector< Real > &gstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
Definition: ROL_HMCR.hpp:238
void checkInputs(void) const
Definition: ROL_HMCR.hpp:112
void getHessVec(Vector< Real > &hv, std::vector< Real > &hvstat, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
Definition: ROL_HMCR.hpp:304
Real lambda_
Definition: ROL_HMCR.hpp:80
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
Definition: ROL_HMCR.hpp:214
Real pnorm_
Definition: ROL_HMCR.hpp:90
Real prob_
Definition: ROL_HMCR.hpp:79
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
Definition: ROL_HMCR.hpp:201
Real coeff0_
Definition: ROL_HMCR.hpp:91
void updateHessVec(Objective< Real > &obj, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for Hessian-time-a-vector computation.
Definition: ROL_HMCR.hpp:266
Real coeff2_
Definition: ROL_HMCR.hpp:93
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
Definition: ROL_HMCR.hpp:188
ROL::Ptr< PlusFunction< Real > > plusFunction_
Definition: ROL_HMCR.hpp:78
Provides the interface to evaluate objective functions.
Provides the interface to implement any functional that maps a random variable to a (extended) real n...
Real computeValue(Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > g_
virtual void initialize(const Vector< Real > &x)
Initialize temporary variables.
void computeHessVec(Vector< Real > &hv, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > hv_
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > dualVector_
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
void sumAll(Real *input, Real *output, int dim) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void scale(const Real alpha)=0
Compute where .
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
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153