ROL
ROL_DouglasRachfordProjection_Def.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
45#ifndef ROL_DOUGLASRACHFORDPROJECTION_DEF_H
46#define ROL_DOUGLASRACHFORDPROJECTION_DEF_H
47
48namespace ROL {
49
50template<typename Real>
52 const Vector<Real> &xdual,
53 const Ptr<BoundConstraint<Real>> &bnd,
54 const Ptr<Constraint<Real>> &con,
55 const Vector<Real> &mul,
56 const Vector<Real> &res)
57 : PolyhedralProjection<Real>(xprim,xdual,bnd,con,mul,res),
58 DEFAULT_atol_ (1e-2*std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
59 DEFAULT_rtol_ (1e-2*std::sqrt(ROL_EPSILON<Real>())),
60 DEFAULT_maxit_ (10000),
61 DEFAULT_verbosity_ (0),
62 DEFAULT_alpha1_ (0.5),
63 DEFAULT_gamma_ (10.0),
64 DEFAULT_t0_ (1.9),
65 atol_ (DEFAULT_atol_),
66 rtol_ (DEFAULT_rtol_),
67 maxit_ (DEFAULT_maxit_),
68 verbosity_ (DEFAULT_verbosity_),
69 alpha1_ (DEFAULT_alpha1_),
70 alpha2_ (1.0-alpha1_),
71 gamma_ (DEFAULT_gamma_),
72 t0_ (DEFAULT_t0_) {
73 dim_ = mul.dimension();
74 tmp_ = xprim.clone();
75 y_ = xprim.clone();
76 q_ = xprim.clone();
77 p_ = xprim.clone();
78 z_ = xdual.clone();
79 if (dim_ == 1) {
80 Real tol(std::sqrt(ROL_EPSILON<Real>()));
81 xprim_->zero();
83 con_->value(*res_,*xprim_,tol);
84 b_ = res_->dot(*res_->basis(0));
85 mul_->setScalar(static_cast<Real>(1));
86 con_->applyAdjointJacobian(*z_,*mul_,xprim,tol);
87 xprim_->set(z_->dual());
88 cdot_ = xprim_->dot(*xprim_);
89 }
90 z_->zero();
91}
92
93template<typename Real>
95 const Vector<Real> &xdual,
96 const Ptr<BoundConstraint<Real>> &bnd,
97 const Ptr<Constraint<Real>> &con,
98 const Vector<Real> &mul,
99 const Vector<Real> &res,
100 ParameterList &list)
101 : DouglasRachfordProjection<Real>(xprim,xdual,bnd,con,mul,res) {
102 atol_ = list.sublist("General").sublist("Polyhedral Projection").get("Absolute Tolerance", DEFAULT_atol_);
103 rtol_ = list.sublist("General").sublist("Polyhedral Projection").get("Relative Tolerance", DEFAULT_rtol_);
104 maxit_ = list.sublist("General").sublist("Polyhedral Projection").get("Iteration Limit", DEFAULT_maxit_);
105 verbosity_ = list.sublist("General").get("Output Level", DEFAULT_verbosity_);
106 alpha1_ = list.sublist("General").sublist("Polyhedral Projection").sublist("Douglas-Rachford").get("Constraint Weight", DEFAULT_alpha1_);
107 alpha2_ = static_cast<Real>(1)-alpha1_;
108 gamma_ = list.sublist("General").sublist("Polyhedral Projection").sublist("Douglas-Rachford").get("Penalty Parameter", DEFAULT_gamma_);
109 t0_ = list.sublist("General").sublist("Polyhedral Projection").sublist("Douglas-Rachford").get("Relaxation Parameter", DEFAULT_t0_);
110}
111
112template<typename Real>
114 if (con_ == nullPtr) {
115 bnd_->project(x);
116 }
117 else {
118 project_DouglasRachford(x, stream);
119 }
120}
121
122template<typename Real>
124 return xprim_->dot(x) + b_;
125}
126
127template<typename Real>
129 Real tol(std::sqrt(ROL_EPSILON<Real>()));
130 con_->update(y,UpdateType::Temp);
131 con_->value(r,y,tol);
132}
133
134template<typename Real>
136 x.set(y);
137 bnd_->project(x);
138}
139
140template<typename Real>
142 if (dim_ == 1) {
143 Real rhs = residual_1d(y);
144 Real lam = -rhs/cdot_;
145 x.set(y);
146 x.axpy(lam,*xprim_);
147 }
148 else {
149 Real tol = std::sqrt(ROL_EPSILON<Real>());
150 residual_nd(*res_,y);
151 con_->solveAugmentedSystem(x,*mul_,*z_,*res_,y,tol);
152 x.scale(static_cast<Real>(-1));
153 x.plus(y);
154 }
155}
156
157template<typename Real>
159 const Real one(1), two(2), xnorm(x.norm()), ctol(std::min(atol_,rtol_*xnorm));
160 Real rnorm(0);
161 p_->zero(); q_->zero(); y_->set(x);
162 std::ios_base::fmtflags streamFlags(stream.flags());
163 if (verbosity_ > 2) {
164 stream << std::scientific << std::setprecision(6);
165 stream << std::endl;
166 stream << " Polyhedral Projection using Douglas Rachford Splitting" << std::endl;
167 stream << " ";
168 stream << std::setw(6) << std::left << "iter";
169 stream << std::setw(15) << std::left << "error";
170 stream << std::setw(15) << std::left << "tol";
171 stream << std::endl;
172 }
173 for (int cnt=0; cnt < maxit_; ++cnt) {
174 // Constraint projection
175 tmp_->set(*y_);
176 tmp_->axpy(alpha1_*gamma_,x);
177 tmp_->scale(one/(alpha1_*gamma_+one));
178 project_con(*p_,*tmp_);
179 // Bounds projection
180 tmp_->zero();
181 tmp_->axpy(two,*p_);
182 tmp_->axpy(-one,*y_);
183 tmp_->axpy(alpha2_*gamma_,x);
184 tmp_->scale(one/(alpha2_*gamma_+one));
185 project_bnd(*q_,*tmp_);
186 // Dual update
187 tmp_->set(*q_);
188 tmp_->axpy(-one,*p_);
189 rnorm = tmp_->norm();
190 if (verbosity_ > 2) {
191 stream << " ";
192 stream << std::setw(6) << std::left << cnt;
193 stream << std::setw(15) << std::left << rnorm;
194 stream << std::setw(15) << std::left << ctol;
195 stream << std::endl;
196 }
197 if (rnorm <= ctol) break;
198 y_->axpy(t0_,*tmp_);
199 }
200 if (verbosity_ > 2) stream << std::endl;
201 x.set(*q_);
202 if (rnorm > ctol) {
203 //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
204 stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
205 stream << rnorm << " rtol = " << ctol << std::endl;
206 }
207 stream.flags(streamFlags);
208}
209
210} // namespace ROL
211
212#endif
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Real residual_1d(const Vector< Real > &x) const
DouglasRachfordProjection(const Vector< Real > &xprim, const Vector< Real > &xdual, const Ptr< BoundConstraint< Real > > &bnd, const Ptr< Constraint< Real > > &con, const Vector< Real > &mul, const Vector< Real > &res)
void project_DouglasRachford(Vector< Real > &x, std::ostream &stream=std::cout) const
void project_bnd(Vector< Real > &x, const Vector< Real > &y) const
void project_con(Vector< Real > &x, const Vector< Real > &y) const
void residual_nd(Vector< Real > &r, const Vector< Real > &y) const
const Ptr< Constraint< Real > > con_
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition ROL_Types.hpp:91