ROL
ROL_DykstraProjection_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_DYKSTRAPROJECTION_DEF_H
46#define ROL_DYKSTRAPROJECTION_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_ (std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
59 DEFAULT_rtol_ (std::sqrt(ROL_EPSILON<Real>())),
60 DEFAULT_maxit_ (10000),
61 DEFAULT_verbosity_ (0),
62 atol_ (DEFAULT_atol_),
63 rtol_ (DEFAULT_rtol_),
64 maxit_ (DEFAULT_maxit_),
65 verbosity_ (DEFAULT_verbosity_) {
66 dim_ = mul.dimension();
67 tmp_ = xprim.clone();
68 y_ = xprim.clone();
69 q_ = xprim.clone();
70 p_ = xprim.clone();
71 z_ = xdual.clone();
72 if (dim_ == 1) {
73 Real tol(std::sqrt(ROL_EPSILON<Real>()));
74 xprim_->zero();
76 con_->value(*res_,*xprim_,tol);
77 b_ = res_->dot(*res_->basis(0));
78 mul_->setScalar(static_cast<Real>(1));
79 con_->applyAdjointJacobian(*z_,*mul_,xprim,tol);
80 xprim_->set(z_->dual());
81 cdot_ = xprim_->dot(*xprim_);
82 }
83 z_->zero();
84}
85
86template<typename Real>
88 const Vector<Real> &xdual,
89 const Ptr<BoundConstraint<Real>> &bnd,
90 const Ptr<Constraint<Real>> &con,
91 const Vector<Real> &mul,
92 const Vector<Real> &res,
93 ParameterList &list)
94 : DykstraProjection<Real>(xprim,xdual,bnd,con,mul,res) {
95 atol_ = list.sublist("General").sublist("Polyhedral Projection").get("Absolute Tolerance", DEFAULT_atol_);
96 rtol_ = list.sublist("General").sublist("Polyhedral Projection").get("Relative Tolerance", DEFAULT_rtol_);
97 maxit_ = list.sublist("General").sublist("Polyhedral Projection").get("Iteration Limit", DEFAULT_maxit_);
98 verbosity_ = list.sublist("General").get("Output Level", DEFAULT_verbosity_);
99}
100
101template<typename Real>
102void DykstraProjection<Real>::project(Vector<Real> &x, std::ostream &stream) {
103 if (con_ == nullPtr) {
104 bnd_->project(x);
105 }
106 else {
107 project_Dykstra(x, stream);
108 }
109}
110
111template<typename Real>
113 return xprim_->dot(x) + b_;
114}
115
116template<typename Real>
118 Real tol(std::sqrt(ROL_EPSILON<Real>()));
119 con_->update(y,UpdateType::Temp);
120 con_->value(r,y,tol);
121}
122
123template<typename Real>
125 x.set(y);
126 bnd_->project(x);
127}
128
129template<typename Real>
131 if (dim_ == 1) {
132 Real rhs = residual_1d(y);
133 Real lam = -rhs/cdot_;
134 x.set(y);
135 x.axpy(lam,*xprim_);
136 }
137 else {
138 Real tol = std::sqrt(ROL_EPSILON<Real>());
139 residual_nd(*res_,y);
140 con_->solveAugmentedSystem(x,*mul_,*z_,*res_,y,tol);
141 x.scale(static_cast<Real>(-1));
142 x.plus(y);
143 }
144}
145
146template<typename Real>
147void DykstraProjection<Real>::project_Dykstra(Vector<Real> &x, std::ostream &stream) const {
148 const Real one(1), xnorm(x.norm()), ctol(std::min(atol_,rtol_*xnorm));
149 Real norm1(0), norm2(0), rnorm(0);
150 p_->zero(); q_->zero();
151 std::ios_base::fmtflags streamFlags(stream.flags());
152 if (verbosity_ > 2) {
153 stream << std::scientific << std::setprecision(6);
154 stream << std::endl;
155 stream << " Polyhedral Projection using Dykstra's Algorithm" << std::endl;
156 stream << " ";
157 stream << std::setw(6) << std::left << "iter";
158 stream << std::setw(15) << std::left << "con norm";
159 stream << std::setw(15) << std::left << "bnd norm";
160 stream << std::setw(15) << std::left << "error";
161 stream << std::setw(15) << std::left << "tol";
162 stream << std::endl;
163 }
164 for (int cnt=0; cnt < maxit_; ++cnt) {
165 // Constraint projection
166 tmp_->set(x); tmp_->plus(*p_);
167 project_con(*y_,*tmp_);
168 p_->set(*tmp_); p_->axpy(-one,*y_);
169 // compute error between pnew and pold
170 tmp_->set(x); tmp_->axpy(-one,*y_);
171 norm1 = tmp_->norm();
172 // Bounds projection
173 tmp_->set(*y_); tmp_->plus(*q_);
174 project_bnd(x,*tmp_);
175 q_->set(*tmp_); q_->axpy(-one,x);
176 // compute error between qnew and qold
177 tmp_->set(x); tmp_->axpy(-one,*y_);
178 norm2 = tmp_->norm();
179 // stopping condition based on Birgin/Raydan paper
180 // Robust Stopping Criteria for Dykstra's Algorithm
181 // SISC Vol. 26, No. 4, 2005
182 rnorm = std::sqrt(norm1*norm1 + norm2*norm2);
183 if (verbosity_ > 2) {
184 stream << " ";
185 stream << std::setw(6) << std::left << cnt;
186 stream << std::setw(15) << std::left << norm1;
187 stream << std::setw(15) << std::left << norm2;
188 stream << std::setw(15) << std::left << rnorm;
189 stream << std::setw(15) << std::left << ctol;
190 stream << std::endl;
191 }
192 if (rnorm <= ctol) break;
193 }
194 if (verbosity_ > 2) {
195 stream << std::endl;
196 }
197 if (rnorm > ctol) {
198 //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
199 stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
200 stream << rnorm << " rtol = " << ctol << std::endl;
201 }
202 stream.flags(streamFlags);
203}
204
205} // namespace ROL
206
207#endif
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
Ptr< Vector< Real > > y_
Ptr< Vector< Real > > z_
Real residual_1d(const Vector< Real > &x) const
void residual_nd(Vector< Real > &r, const Vector< Real > &y) const
Ptr< Vector< Real > > p_
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
void project_Dykstra(Vector< Real > &x, std::ostream &stream=std::cout) const
Ptr< Vector< Real > > tmp_
void project_con(Vector< Real > &x, const Vector< Real > &y) const
DykstraProjection(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_bnd(Vector< Real > &x, const Vector< Real > &y) const
Ptr< Vector< Real > > q_
const Ptr< Constraint< Real > > con_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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.
Definition: ROL_Vector.hpp:196
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91