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
25namespace ROL {
26
27template <class Real>
28class PrimalDualRisk {
29private:
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
68public:
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
247private:
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_PD_BPOE.hpp
ROL_PD_CVaR.hpp
ROL_PD_HMCR2.hpp
ROL_PD_MeanSemiDeviationFromTarget.hpp
ROL_PD_MeanSemiDeviation.hpp
ROL_RiskBoundConstraint.hpp
ROL_RiskLessConstraint.hpp
ROL_RiskVector.hpp
ROL_Solver.hpp
ROL_StochasticObjective.hpp
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::PrimalDualRisk::pd_objective_
Ptr< StochasticObjective< Real > > pd_objective_
Definition ROL_PrimalDualRisk.hpp:56
ROL::PrimalDualRisk::input_
const Ptr< Problem< Real > > input_
Definition ROL_PrimalDualRisk.hpp:30
ROL::PrimalDualRisk::pd_linear_constraint_
Ptr< Constraint< Real > > pd_linear_constraint_
Definition ROL_PrimalDualRisk.hpp:60
ROL::PrimalDualRisk::checkStatus
bool checkStatus(const AlgorithmState< Real > &state, std::ostream &outStream) const
Definition ROL_PrimalDualRisk.hpp:305
ROL::PrimalDualRisk::check
void check(std::ostream &outStream=std::cout)
Definition ROL_PrimalDualRisk.hpp:187
ROL::PrimalDualRisk::pd_vector_
Ptr< Vector< Real > > pd_vector_
Definition ROL_PrimalDualRisk.hpp:57
ROL::PrimalDualRisk::pd_problem_
Ptr< Problem< Real > > pd_problem_
Definition ROL_PrimalDualRisk.hpp:61
ROL::PrimalDualRisk::PrimalDualRisk
PrimalDualRisk(const Ptr< Problem< Real > > &input, const Ptr< SampleGenerator< Real > > &sampler, ParameterList &parlist)
Definition ROL_PrimalDualRisk.hpp:69
ROL::PrimalDualRisk::rvf_
Ptr< PD_RandVarFunctional< Real > > rvf_
Definition ROL_PrimalDualRisk.hpp:32
ROL::PrimalDualRisk::printHeader
void printHeader(std::ostream &outStream) const
Definition ROL_PrimalDualRisk.hpp:248
ROL::PrimalDualRisk::pd_constraint_
Ptr< Constraint< Real > > pd_constraint_
Definition ROL_PrimalDualRisk.hpp:59
ROL::PrimalDualRisk::sampler_
const Ptr< SampleGenerator< Real > > sampler_
Definition ROL_PrimalDualRisk.hpp:31
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::PrimalDualRisk::pd_bound_
Ptr< BoundConstraint< Real > > pd_bound_
Definition ROL_PrimalDualRisk.hpp:58
ROL::RiskVector::getVector
Ptr< const Vector< Real > > getVector(void) const
Definition ROL_RiskVector.hpp:463
ROL::ROL_INF
Real ROL_INF(void)
Definition ROL_Types.hpp:71
ROL::EXITSTATUS_CONVERGED
@ EXITSTATUS_CONVERGED
Definition ROL_Types.hpp:84
ROL::EXITSTATUS_USERDEFINED
@ EXITSTATUS_USERDEFINED
Definition ROL_Types.hpp:88
ROL::AlgorithmState
State for algorithm class. Will be used for restarts.
Definition ROL_Types.hpp:109
Generated by doxygen 1.9.8