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