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"
51#include "ROL_PD_CVaR.hpp"
52#include "ROL_PD_BPOE.hpp"
53#include "ROL_PD_HMCR2.hpp"
54#include "ROL_RiskVector.hpp"
57#include "ROL_Types.hpp"
58
59namespace ROL {
60
61template <class Real>
63private:
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_;
77 Real lalpha_;
78 Real lbeta_;
79 // Subproblem solver tolerances
80 Real gtol_;
81 Real ctol_;
82 Real ltol_;
83 // Penalty parameter information
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
99 Real lnorm_;
100 std::string name_;
101
102public:
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_);
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)) {
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
281private:
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
Contains definitions of custom data types in ROL.
Ptr< StochasticObjective< Real > > pd_objective_
const Ptr< Problem< Real > > input_
Ptr< Constraint< Real > > pd_linear_constraint_
bool checkStatus(const AlgorithmState< Real > &state, std::ostream &outStream) const
void check(std::ostream &outStream=std::cout)
Ptr< Vector< Real > > pd_vector_
Ptr< Problem< Real > > pd_problem_
PrimalDualRisk(const Ptr< Problem< Real > > &input, const Ptr< SampleGenerator< Real > > &sampler, ParameterList &parlist)
Ptr< PD_RandVarFunctional< Real > > rvf_
void printHeader(std::ostream &outStream) const
Ptr< Constraint< Real > > pd_constraint_
const Ptr< SampleGenerator< Real > > sampler_
void run(std::ostream &outStream=std::cout)
void print(const AlgorithmState< Real > &state, std::ostream &outStream) const
Ptr< BoundConstraint< Real > > pd_bound_
Ptr< const Vector< Real > > getVector(void) const
Real ROL_INF(void)
Definition: ROL_Types.hpp:105
@ EXITSTATUS_CONVERGED
Definition: ROL_Types.hpp:118
@ EXITSTATUS_USERDEFINED
Definition: ROL_Types.hpp:122
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:143