ROL
ROL_PrimalDualRisk.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_PRIMALDUALRISK_H
45 #define ROL_PRIMALDUALRISK_H
46 
47 #include "ROL_Solver.hpp"
48 #include "ROL_StochasticObjective.hpp"
49 #include "ROL_PD_MeanSemiDeviation.hpp"
50 #include "ROL_PD_MeanSemiDeviationFromTarget.hpp"
51 #include "ROL_PD_CVaR.hpp"
52 #include "ROL_PD_BPOE.hpp"
53 #include "ROL_PD_HMCR2.hpp"
54 #include "ROL_RiskVector.hpp"
55 #include "ROL_RiskBoundConstraint.hpp"
56 #include "ROL_RiskLessConstraint.hpp"
57 #include "ROL_Types.hpp"
58 
59 namespace ROL {
60 
61 template <class Real>
62 class PrimalDualRisk {
63 private:
64  const Ptr<Problem<Real>> input_;
65  const Ptr<SampleGenerator<Real>> sampler_;
66  Ptr<PD_RandVarFunctional<Real>> rvf_;
67  ParameterList parlist_;
68  // General algorithmic parameters
69  int maxit_;
70  bool print_;
71  Real gtolmin_;
72  Real ctolmin_;
73  Real ltolmin_;
74  Real ltolupdate_;
75  Real tolupdate0_;
76  Real tolupdate1_;
77  Real lalpha_;
78  Real lbeta_;
79  // Subproblem solver tolerances
80  Real gtol_;
81  Real ctol_;
82  Real ltol_;
83  // Penalty parameter information
84  Real penaltyParam_;
85  Real maxPen_;
86  // Forced udpate information
87  Real update_;
88  int freq_;
89 
90  Ptr<StochasticObjective<Real>> pd_objective_;
91  Ptr<Vector<Real>> pd_vector_;
92  Ptr<BoundConstraint<Real>> pd_bound_;
93  Ptr<Constraint<Real>> pd_constraint_;
94  Ptr<Constraint<Real>> pd_linear_constraint_;
95  Ptr<Problem<Real>> pd_problem_;
96 
97  int iter_, nfval_, ngrad_, ncval_;
98  bool converged_;
99  Real lnorm_;
100  std::string name_;
101 
102 public:
103  PrimalDualRisk(const Ptr<Problem<Real>> &input,
104  const Ptr<SampleGenerator<Real>> &sampler,
105  ParameterList &parlist)
106  : input_(input), sampler_(sampler), parlist_(parlist),
107  iter_(0), converged_(true), lnorm_(ROL_INF<Real>()) {
108  // Get status test information
109  maxit_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Iteration Limit",100);
110  print_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Print Subproblem Solve History",false);
111  gtolmin_ = parlist.sublist("Status Test").get("Gradient Tolerance", 1e-8);
112  ctolmin_ = parlist.sublist("Status Test").get("Constraint Tolerance", 1e-8);
113  ltolmin_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Dual Tolerance",1e-6);
114  gtolmin_ = (gtolmin_ <= static_cast<Real>(0) ? std::sqrt(ROL_EPSILON<Real>()) : gtolmin_);
115  ctolmin_ = (ctolmin_ <= static_cast<Real>(0) ? std::sqrt(ROL_EPSILON<Real>()) : ctolmin_);
116  ltolmin_ = (ltolmin_ <= static_cast<Real>(0) ? 1e2*std::sqrt(ROL_EPSILON<Real>()) : ltolmin_);
117  // Get solver tolerances
118  gtol_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Initial Gradient Tolerance", 1e-4);
119  ctol_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Initial Constraint Tolerance", 1e-4);
120  ltol_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Initial Dual Tolerance", 1e-2);
121  ltolupdate_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Dual Tolerance Update Scale", 1e-1);
122  tolupdate0_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Solver Tolerance Decrease Scale", 9e-1);
123  tolupdate1_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Solver Tolerance Update Scale", 1e-1);
124  lalpha_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Dual Tolerance Decrease Exponent", 1e-1);
125  lbeta_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Dual Tolerance Update Exponent", 9e-1);
126  gtol_ = std::max(gtol_, gtolmin_);
127  ctol_ = std::max(ctol_, ctolmin_);
128  ltol_ = std::max(ltol_, ltolmin_);
129  // Get penalty parameter
130  penaltyParam_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Initial Penalty Parameter", 10.0);
131  maxPen_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Maximum Penalty Parameter", -1.0);
132  update_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Penalty Update Scale", 10.0);
133  maxPen_ = (maxPen_ <= static_cast<Real>(0) ? ROL_INF<Real>() : maxPen_);
134  penaltyParam_ = std::min(penaltyParam_,maxPen_);
135  // Get update parameters
136  freq_ = parlist.sublist("SOL").sublist("Primal Dual Risk").get("Update Frequency", 0);
137  // Create risk vector and risk-averse objective
138  ParameterList olist; olist.sublist("SOL") = parlist.sublist("SOL").sublist("Objective");
139  std::string type = olist.sublist("SOL").get<std::string>("Type");
140  if (type == "Risk Averse") {
141  name_ = olist.sublist("SOL").sublist("Risk Measure").get<std::string>("Name");
142  }
143  else if (type == "Probability") {
144  name_ = olist.sublist("SOL").sublist("Probability"). get<std::string>("Name");
145  }
146  else {
147  std::stringstream message;
148  message << ">>> " << type << " is not implemented!";
149  throw Exception::NotImplemented(message.str());
150  }
151  Ptr<ParameterList> parlistptr = makePtrFromRef<ParameterList>(olist);
152  if (name_ == "CVaR") {
153  parlistptr->sublist("SOL").set("Type","Risk Averse");
154  parlistptr->sublist("SOL").sublist("Risk Measure").set("Name","CVaR");
155  Real alpha = olist.sublist("SOL").sublist("Risk Measure").sublist("CVaR").get("Convex Combination Parameter", 1.0);
156  Real beta = olist.sublist("SOL").sublist("Risk Measure").sublist("CVaR").get("Confidence Level", 0.9);
157  rvf_ = makePtr<PD_CVaR<Real>>(alpha, beta);
158  }
159  else if (name_ == "Mean Plus Semi-Deviation") {
160  parlistptr->sublist("SOL").set("Type","Risk Averse");
161  parlistptr->sublist("SOL").sublist("Risk Measure").set("Name","Mean Plus Semi-Deviation");
162  Real coeff = olist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Semi-Deviation").get("Coefficient", 0.5);
163  rvf_ = makePtr<PD_MeanSemiDeviation<Real>>(coeff);
164  }
165  else if (name_ == "Mean Plus Semi-Deviation From Target") {
166  parlistptr->sublist("SOL").set("Type","Risk Averse");
167  parlistptr->sublist("SOL").sublist("Risk Measure").set("Name","Mean Plus Semi-Deviation From Target");
168  Real coeff = olist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Semi-Deviation From Target").get("Coefficient", 0.5);
169  Real target = olist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Semi-Deviation From Target").get("Target", 1.0);
170  rvf_ = makePtr<PD_MeanSemiDeviationFromTarget<Real>>(coeff, target);
171  }
172  else if (name_ == "HMCR") {
173  parlistptr->sublist("SOL").set("Type","Risk Averse");
174  parlistptr->sublist("SOL").sublist("Risk Measure").set("Name","HMCR");
175  //Real alpha = olist.sublist("SOL").sublist("Risk Measure").sublist("HMCR").get("Convex Combination Parameter", 1.0);
176  Real beta = olist.sublist("SOL").sublist("Risk Measure").sublist("HMCR").get("Confidence Level", 0.9);
177  rvf_ = makePtr<PD_HMCR2<Real>>(beta);
178  }
179  else if (name_ == "bPOE") {
180  parlistptr->sublist("SOL").set("Type","Probability");
181  parlistptr->sublist("SOL").sublist("Probability").set("Name","bPOE");
182  Real thresh = olist.sublist("SOL").sublist("Probability").sublist("bPOE").get("Threshold", 1.0);
183  rvf_ = makePtr<PD_BPOE<Real>>(thresh);
184  }
185  else {
186  std::stringstream message;
187  message << ">>> " << name_ << " is not implemented!";
188  throw Exception::NotImplemented(message.str());
189  }
190  pd_vector_ = makePtr<RiskVector<Real>>(parlistptr,
191  input_->getPrimalOptimizationVector());
192  rvf_->setData(*sampler_, penaltyParam_);
193  pd_objective_ = makePtr<StochasticObjective<Real>>(input_->getObjective(),
194  rvf_, sampler_, true);
195  // Create risk bound constraint
196  pd_bound_ = makePtr<RiskBoundConstraint<Real>>(parlistptr,
197  input_->getBoundConstraint());
198  // Create riskless constraint
199  pd_constraint_ = nullPtr;
200  if (input_->getConstraint() != nullPtr) {
201  pd_constraint_ = makePtr<RiskLessConstraint<Real>>(input_->getConstraint());
202  }
203  pd_linear_constraint_ = nullPtr;
204  if (input_->getPolyhedralProjection() != nullPtr) {
205  pd_linear_constraint_ = makePtr<RiskLessConstraint<Real>>(input_->getPolyhedralProjection()->getLinearConstraint());
206  }
207  // Build primal-dual subproblems
208  pd_problem_ = makePtr<Problem<Real>>(pd_objective_, pd_vector_);
209  if (pd_bound_->isActivated()) {
210  pd_problem_->addBoundConstraint(pd_bound_);
211  }
212  if (pd_constraint_ != nullPtr) {
213  pd_problem_->addConstraint("PD Constraint",pd_constraint_,input_->getMultiplierVector());
214  }
215  if (pd_linear_constraint_ != nullPtr) {
216  pd_problem_->addLinearConstraint("PD Linear Constraint",pd_linear_constraint_,input_->getPolyhedralProjection()->getMultiplier());
217  pd_problem_->setProjectionAlgorithm(parlist);
218  }
219  }
220 
221  void check(std::ostream &outStream = std::cout) {
222  pd_problem_->check(true,outStream);
223  }
224 
225  void run(std::ostream &outStream = std::cout) {
226  const Real one(1);
227  Real theta(1);
228  int spiter(0);
229  iter_ = 0; converged_ = true; lnorm_ = ROL_INF<Real>();
230  nfval_ = 0; ncval_ = 0; ngrad_ = 0;
231  // Output
232  printHeader(outStream);
233  Ptr<Solver<Real>> solver;
234  for (iter_ = 0; iter_ < maxit_; ++iter_) {
235  parlist_.sublist("Status Test").set("Gradient Tolerance", gtol_);
236  parlist_.sublist("Status Test").set("Constraint Tolerance", ctol_);
237  solver = makePtr<Solver<Real>>(pd_problem_, parlist_);
238  if (print_) solver->solve(outStream);
239  else solver->solve();
240  converged_ = (solver->getAlgorithmState()->statusFlag == EXITSTATUS_CONVERGED
241  ||solver->getAlgorithmState()->statusFlag == EXITSTATUS_USERDEFINED
242  ? true : false);
243  spiter += solver->getAlgorithmState()->iter;
244  nfval_ += solver->getAlgorithmState()->nfval;
245  ngrad_ += solver->getAlgorithmState()->ngrad;
246  ncval_ += solver->getAlgorithmState()->ncval;
247  lnorm_ = rvf_->computeDual(*sampler_);
248  // Output
249  print(*solver->getAlgorithmState(),outStream);
250  // Check termination conditions
251  if (checkStatus(*solver->getAlgorithmState(),outStream)) break;
252  // Update penalty parameter and solver tolerances
253  rvf_->updateDual(*sampler_);
254  if (converged_) {
255  if (lnorm_ > penaltyParam_*ltol_ || (freq_ > 0 && iter_%freq_ == 0)) {
256  penaltyParam_ = std::min(update_*penaltyParam_, maxPen_);
257  rvf_->updatePenalty(penaltyParam_);
258  theta = std::min(one/penaltyParam_,one);
259  ltol_ = std::max(ltolupdate_*std::pow(theta,lalpha_), ltolmin_);
260  gtol_ = std::max(tolupdate0_*gtol_, gtolmin_);
261  ctol_ = std::max(tolupdate0_*ctol_, ctolmin_);
262  }
263  else {
264  theta = std::min(one/penaltyParam_,one);
265  ltol_ = std::max(ltol_*std::pow(theta,lbeta_), ltolmin_);
266  gtol_ = std::max(tolupdate1_*gtol_, gtolmin_);
267  ctol_ = std::max(tolupdate1_*ctol_, ctolmin_);
268  }
269  }
270  }
271  input_->getPrimalOptimizationVector()->set(
272  *dynamicPtrCast<RiskVector<Real>>(pd_problem_->getPrimalOptimizationVector())->getVector());
273  // Output reason for termination
274  if (iter_ >= maxit_) {
275  outStream << "Maximum number of iterations exceeded" << std::endl;
276  }
277  outStream << "Primal Dual Risk required " << spiter
278  << " subproblem iterations" << std::endl;
279  }
280 
281 private:
282  void printHeader(std::ostream &outStream) const {
283  std::ios_base::fmtflags flags = outStream.flags();
284  outStream << std::scientific << std::setprecision(6);
285  outStream << std::endl << "Primal Dual Risk Minimization using "
286  << name_ << std::endl << " "
287  << std::setw(8) << std::left << "iter"
288  << std::setw(15) << std::left << "value";
289  if (pd_constraint_ != nullPtr) {
290  outStream << std::setw(15) << std::left << "cnorm";
291  }
292  outStream << std::setw(15) << std::left << "gnorm"
293  << std::setw(15) << std::left << "lnorm"
294  << std::setw(15) << std::left << "penalty";
295  if (pd_constraint_ != nullPtr) {
296  outStream << std::setw(15) << std::left << "ctol";
297  }
298  outStream << std::setw(15) << std::left << "gtol"
299  << std::setw(15) << std::left << "ltol"
300  << std::setw(10) << std::left << "nfval"
301  << std::setw(10) << std::left << "ngrad";
302  if (pd_constraint_ != nullPtr) {
303  outStream << std::setw(10) << std::left << "ncval";
304  }
305  outStream << std::setw(10) << std::left << "subiter"
306  << std::setw(8) << std::left << "success"
307  << std::endl;
308  outStream.setf(flags);
309  }
310 
311  void print(const AlgorithmState<Real> &state, std::ostream &outStream) const {
312  std::ios_base::fmtflags flags = outStream.flags();
313  outStream << std::scientific << std::setprecision(6);
314  outStream << " "
315  << std::setw(8) << std::left << iter_+1
316  << std::setw(15) << std::left << state.value;
317  if (pd_constraint_ != nullPtr) {
318  outStream << std::setw(15) << std::left << state.cnorm;
319  }
320  outStream << std::setw(15) << std::left << state.gnorm
321  << std::setw(15) << std::left << lnorm_
322  << std::setw(15) << std::left << penaltyParam_;
323  if (pd_constraint_ != nullPtr) {
324  outStream << std::setw(15) << std::left << ctol_;
325  }
326  outStream << std::setw(15) << std::left << gtol_
327  << std::setw(15) << std::left << ltol_
328  << std::setw(10) << std::left << nfval_
329  << std::setw(10) << std::left << ngrad_;
330  if (pd_constraint_ != nullPtr) {
331  outStream << std::setw(10) << std::left << ncval_;
332  }
333  outStream << std::setw(10) << std::left << state.iter
334  << std::setw(8) << std::left << converged_
335  << std::endl;
336  outStream.setf(flags);
337  }
338 
339  bool checkStatus(const AlgorithmState<Real> &state, std::ostream &outStream) const {
340  bool flag = false;
341 // if (converged_ && state.iter==0 && lnorm_ < tol) {
342 // outStream << "Subproblem solve converged in zero iterations"
343 // << " and the difference in the multipliers was less than "
344 // << tol1 << std::endl;
345 // flag = true;
346 // }
347  if (pd_constraint_ == nullPtr) {
348  if (state.gnorm < gtolmin_ && lnorm_/penaltyParam_ < ltolmin_) {
349  outStream << "Solver tolerance met"
350  << " and the difference in the multipliers was less than "
351  << ltolmin_ << std::endl;
352  flag = true;
353  }
354  }
355  else {
356  if (state.gnorm < gtolmin_ && state.cnorm < ctolmin_ && lnorm_/penaltyParam_ < ltolmin_) {
357  outStream << "Solver tolerance met"
358  << " and the difference in the multipliers was less than "
359  << ltolmin_ << std::endl;
360  flag = true;
361  }
362  }
363  return flag;
364  }
365 
366 }; // class PrimalDualRisk
367 
368 } // namespace ROL
369 
370 #endif
ROL::PrimalDualRisk::rvf_
Ptr< PD_RandVarFunctional< Real > > rvf_
Definition: ROL_PrimalDualRisk.hpp:66
ROL_PD_MeanSemiDeviation.hpp
ROL::PrimalDualRisk::input_
const Ptr< Problem< Real > > input_
Definition: ROL_PrimalDualRisk.hpp:64
ROL_StochasticObjective.hpp
ROL::PrimalDualRisk::pd_vector_
Ptr< Vector< Real > > pd_vector_
Definition: ROL_PrimalDualRisk.hpp:91
ROL_RiskLessConstraint.hpp
ROL::PrimalDualRisk::pd_linear_constraint_
Ptr< Constraint< Real > > pd_linear_constraint_
Definition: ROL_PrimalDualRisk.hpp:94
ROL_PD_MeanSemiDeviationFromTarget.hpp
ROL::ROL_INF
Real ROL_INF(void)
Definition: ROL_Types.hpp:105
ROL::PrimalDualRisk::printHeader
void printHeader(std::ostream &outStream) const
Definition: ROL_PrimalDualRisk.hpp:282
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL_PD_CVaR.hpp
ROL::PrimalDualRisk::run
void run(std::ostream &outStream=std::cout)
Definition: ROL_PrimalDualRisk.hpp:225
ROL::PrimalDualRisk::print
void print(const AlgorithmState< Real > &state, std::ostream &outStream) const
Definition: ROL_PrimalDualRisk.hpp:311
ROL_PD_BPOE.hpp
ROL::PrimalDualRisk::PrimalDualRisk
PrimalDualRisk(const Ptr< Problem< Real >> &input, const Ptr< SampleGenerator< Real >> &sampler, ParameterList &parlist)
Definition: ROL_PrimalDualRisk.hpp:103
ROL::AlgorithmState
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:143
ROL::PrimalDualRisk::checkStatus
bool checkStatus(const AlgorithmState< Real > &state, std::ostream &outStream) const
Definition: ROL_PrimalDualRisk.hpp:339
ROL::PrimalDualRisk::pd_problem_
Ptr< Problem< Real > > pd_problem_
Definition: ROL_PrimalDualRisk.hpp:95
ROL_RiskVector.hpp
ROL::PrimalDualRisk::pd_constraint_
Ptr< Constraint< Real > > pd_constraint_
Definition: ROL_PrimalDualRisk.hpp:93
ROL_PD_HMCR2.hpp
ROL::PrimalDualRisk::sampler_
const Ptr< SampleGenerator< Real > > sampler_
Definition: ROL_PrimalDualRisk.hpp:65
ROL::PrimalDualRisk::pd_objective_
Ptr< StochasticObjective< Real > > pd_objective_
Definition: ROL_PrimalDualRisk.hpp:90
ROL::PrimalDualRisk::check
void check(std::ostream &outStream=std::cout)
Definition: ROL_PrimalDualRisk.hpp:221
ROL::PrimalDualRisk::pd_bound_
Ptr< BoundConstraint< Real > > pd_bound_
Definition: ROL_PrimalDualRisk.hpp:92
ROL_RiskBoundConstraint.hpp
ROL_Solver.hpp
Generated by   doxygen 1.8.5