ROL
ROL_Bundle_U_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#ifndef ROL_BUNDLE_U_DEF_H
45#define ROL_BUNDLE_U_DEF_H
46
47#include "ROL_Types.hpp"
48
49namespace ROL {
50
51template<typename Real>
52void Bundle_U<Real>::remove(const std::vector<unsigned> &ind) {
53 Real zero(0);
54 for (unsigned j = ind.back()+1; j < size_; ++j) {
55 (subgradients_[j-1])->set(*(subgradients_[j]));
56 linearizationErrors_[j-1] = linearizationErrors_[j];
57 distanceMeasures_[j-1] = distanceMeasures_[j];
58 dualVariables_[j-1] = dualVariables_[j];
59 }
60 (subgradients_[size_-1])->zero();
61 linearizationErrors_[size_-1] = ROL_OVERFLOW<Real>();
62 distanceMeasures_[size_-1] = ROL_OVERFLOW<Real>();
63 dualVariables_[size_-1] = zero;
64 for (unsigned i = ind.size()-1; i > 0; --i) {
65 for (unsigned j = ind[i-1]+1; j < size_; ++j) {
66 (subgradients_[j-1])->set(*(subgradients_[j]));
67 linearizationErrors_[j-1] = linearizationErrors_[j];
68 distanceMeasures_[j-1] = distanceMeasures_[j];
69 dualVariables_[j-1] = dualVariables_[j];
70 }
71 }
72 size_ -= ind.size();
73}
74
75template<typename Real>
76void Bundle_U<Real>::add(const Vector<Real> &g, const Real le, const Real dm) {
77 Real zero(0);
78 (subgradients_[size_])->set(g);
79 linearizationErrors_[size_] = le;
80 distanceMeasures_[size_] = dm;
81 dualVariables_[size_] = zero;
82 size_++;
83}
84
85template<typename Real>
86Bundle_U<Real>::Bundle_U(const unsigned maxSize,
87 const Real coeff,
88 const Real omega,
89 const unsigned remSize)
90 : size_(0), maxSize_(maxSize), remSize_(remSize),
91 coeff_(coeff), omega_(omega), isInitialized_(false) {
92 Real zero(0);
93 remSize_ = ((remSize_ < 2) ? 2 : ((remSize_ > maxSize_-1) ? maxSize_-1 : remSize_));
94 coeff_ = std::max(static_cast<Real>(0),coeff_);
95 omega_ = std::max(static_cast<Real>(1),omega_);
96 subgradients_.clear();
97 subgradients_.assign(maxSize,nullPtr);
99 linearizationErrors_.assign(maxSize_,ROL_OVERFLOW<Real>());
100 distanceMeasures_.clear();
101 distanceMeasures_.assign(maxSize_,ROL_OVERFLOW<Real>());
102 dualVariables_.clear();
104}
105
106template<typename Real>
108 if ( !isInitialized_ ) {
109 Real zero(0), one(1);
110 for (unsigned i = 0; i < maxSize_; ++i) {
111 subgradients_[i] = g.clone();
112 }
113 (subgradients_[0])->set(g);
114 linearizationErrors_[0] = zero;
115 distanceMeasures_[0] = zero;
116 dualVariables_[0] = one;
117 size_++;
118 isInitialized_ = true;
119 tG_ = g.clone();
120 yG_ = g.clone();
121 eG_ = g.clone();
122 gx_ = g.clone();
123 ge_ = g.clone();
124 }
125}
126
127template<typename Real>
128const Real Bundle_U<Real>::linearizationError(const unsigned i) const {
129 return linearizationErrors_[i];
130}
131
132template<typename Real>
133const Real Bundle_U<Real>::distanceMeasure(const unsigned i) const {
134 return distanceMeasures_[i];
135}
136
137template<typename Real>
138const Vector<Real> & Bundle_U<Real>::subgradient(const unsigned i) const {
139 return *(subgradients_[i]);
140}
141
142template<typename Real>
143const Real Bundle_U<Real>::getDualVariable(const unsigned i) const {
144 return dualVariables_[i];
145}
146
147template<typename Real>
148void Bundle_U<Real>::setDualVariable(const unsigned i, const Real val) {
149 dualVariables_[i] = val;
150}
151
152template<typename Real>
154 const Real zero(0);
155 dualVariables_.assign(size_,zero);
156}
157
158template<typename Real>
159const Real Bundle_U<Real>::computeAlpha(const Real dm, const Real le) const {
160 Real alpha = le;
161 if ( coeff_ > ROL_EPSILON<Real>() ) {
162 alpha = std::max(coeff_*std::pow(dm,omega_),le);
163 }
164 return alpha;
165}
166
167template<typename Real>
168const Real Bundle_U<Real>::alpha(const unsigned i) const {
169 return computeAlpha(distanceMeasures_[i],linearizationErrors_[i]);
170}
171
172template<typename Real>
173unsigned Bundle_U<Real>::size(void) const {
174 return size_;
175}
176
177template<typename Real>
178void Bundle_U<Real>::aggregate(Vector<Real> &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const {
179 Real zero(0), one(1);
180 aggSubGrad.zero(); aggLinErr = zero; aggDistMeas = zero; eG_->zero();
181 Real eLE(0), eDM(0), yLE(0), yDM(0), tLE(0), tDM(0);
182 for (unsigned i = 0; i < size_; ++i) {
183 // Compute aggregate subgradient using Kahan's compensated sum
184 //aggSubGrad.axpy(dualVariables_[i],*subgradients_[i]);
185 yG_->set(*subgradients_[i]); yG_->scale(dualVariables_[i]); yG_->axpy(-one,*eG_);
186 tG_->set(aggSubGrad); tG_->plus(*yG_);
187 eG_->set(*tG_); eG_->axpy(-one,aggSubGrad); eG_->axpy(-one,*yG_);
188 aggSubGrad.set(*tG_);
189 // Compute aggregate linearization error using Kahan's compensated sum
190 //aggLinErr += dualVariables_[i]*linearizationErrors_[i];
191 yLE = dualVariables_[i]*linearizationErrors_[i] - eLE;
192 tLE = aggLinErr + yLE;
193 eLE = (tLE - aggLinErr) - yLE;
194 aggLinErr = tLE;
195 // Compute aggregate distance measure using Kahan's compensated sum
196 //aggDistMeas += dualVariables_[i]*distanceMeasures_[i];
197 yDM = dualVariables_[i]*distanceMeasures_[i] - eDM;
198 tDM = aggDistMeas + yDM;
199 eDM = (tDM - aggDistMeas) - yDM;
200 aggDistMeas = tDM;
201 }
202}
203
204template<typename Real>
205void Bundle_U<Real>::reset(const Vector<Real> &g, const Real le, const Real dm) {
206 if (size_ == maxSize_) {
207 // Find indices to remove
208 unsigned loc = size_, cnt = 0;
209 std::vector<unsigned> ind(remSize_,0);
210 for (unsigned i = size_; i > 0; --i) {
211 if ( std::abs(linearizationErrors_[i-1]) < ROL_EPSILON<Real>() ) {
212 loc = i-1;
213 break;
214 }
215 }
216 for (unsigned i = 0; i < size_; ++i) {
217 if ( i != loc ) {
218 ind[cnt] = i;
219 cnt++;
220 }
221 if (cnt == remSize_) {
222 break;
223 }
224 }
225 // Remove indices
226 remove(ind);
227 // Add aggregate subgradient
228 add(g,le,dm);
229 }
230}
231
232template<typename Real>
233void Bundle_U<Real>::update(const bool flag, const Real linErr, const Real distMeas,
234 const Vector<Real> &g, const Vector<Real> &s) {
235 Real zero(0);
236 if ( flag ) {
237 // Serious step taken: Update linearlization errors and distance measures
238 for (unsigned i = 0; i < size_; ++i) {
239 //linearizationErrors_[i] += linErr - subgradients_[i]->dot(s.dual());
240 linearizationErrors_[i] += linErr - subgradients_[i]->apply(s);
241 distanceMeasures_[i] += distMeas;
242 }
243 linearizationErrors_[size_] = zero;
244 distanceMeasures_[size_] = zero;
245 }
246 else {
247 // Null step taken:
248 linearizationErrors_[size_] = linErr;
249 distanceMeasures_[size_] = distMeas;
250 }
251 // Update (sub)gradient bundle
252 (subgradients_[size_])->set(g);
253 // Update dual variables
254 dualVariables_[size_] = zero;
255 // Update bundle size
256 size_++;
257}
258
259template<typename Real>
260const Real Bundle_U<Real>::GiGj(const unsigned i, const unsigned j) const {
261 return subgradient(i).dot(subgradient(j));
262}
263
264template<typename Real>
265const Real Bundle_U<Real>::dotGi(const unsigned i, const Vector<Real> &x) const {
266 return x.dot(subgradient(i));
267}
268
269template<typename Real>
270void Bundle_U<Real>::addGi(const unsigned i, const Real a, Vector<Real> &x) const {
271 x.axpy(a,subgradient(i));
272}
273
274template<typename Real>
275Real Bundle_U<Real>::evaluateObjective(std::vector<Real> &g, const std::vector<Real> &x, const Real t) const {
276 Real one(1), half(0.5);
277 gx_->zero(); eG_->zero();
278 for (unsigned i = 0; i < size(); ++i) {
279 // Compute Gx using Kahan's compensated sum
280 //gx_->axpy(x[i],*subgradients_[i]);
281 yG_->set(subgradient(i)); yG_->scale(x[i]); yG_->axpy(-one,*eG_);
282 tG_->set(*gx_); tG_->plus(*yG_);
283 eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->axpy(-one,*yG_);
284 gx_->set(*tG_);
285 }
286 Real Hx(0), val(0), err(0), tmp(0), y(0);
287 for (unsigned i = 0; i < size(); ++i) {
288 // Compute < g_i, Gx >
289 Hx = dotGi(i,*gx_);
290 // Add to the objective function value using Kahan's compensated sum
291 //val += x[i]*(half*Hx + alpha(i)/t);
292 y = x[i]*(half*Hx + alpha(i)/t) - err;
293 tmp = val + y;
294 err = (tmp - val) - y;
295 val = tmp;
296 // Add gradient component
297 g[i] = Hx + alpha(i)/t;
298 }
299 return val;
300}
301
302template<typename Real>
303unsigned Bundle_U<Real>::solveDual_dim1(const Real t, const unsigned maxit, const Real tol) {
304 setDualVariable(0,static_cast<Real>(1));
305 //std::cout << "dim = " << Bundle<Real>::size() << " iter = " << 0 << " CONVERGED!\n";
306 return 0;
307}
308
309template<typename Real>
310unsigned Bundle_U<Real>::solveDual_dim2(const Real t, const unsigned maxit, const Real tol) {
311 Real diffg = gx_->dot(*gx_), zero(0), one(1), half(0.5);
312 gx_->set(subgradient(0)); addGi(1,-one,*gx_);
313 if ( std::abs(diffg) > ROL_EPSILON<Real>() ) {
314 Real diffa = (alpha(0)-alpha(1))/t;
315 Real gdiffg = dotGi(1,*gx_);
316 setDualVariable(0,std::min(one,std::max(zero,-(gdiffg+diffa)/diffg)));
317 setDualVariable(1,one-getDualVariable(0));
318 }
319 else {
320 if ( std::abs(alpha(0)-alpha(1)) > ROL_EPSILON<Real>() ) {
321 if ( alpha(0) < alpha(1) ) {
322 setDualVariable(0,one); setDualVariable(1,zero);
323 }
324 else if ( alpha(0) > alpha(1) ) {
325 setDualVariable(0,zero); setDualVariable(1,one);
326 }
327 }
328 else {
329 setDualVariable(0,half); setDualVariable(1,half);
330 }
331 }
332 //std::cout << "dim = " << Bundle<Real>::size() << " iter = " << 0 << " CONVERGED!\n";
333 return 0;
334}
335
336} // namespace ROL
337
338#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of custom data types in ROL.
void update(const bool flag, const Real linErr, const Real distMeas, const Vector< Real > &g, const Vector< Real > &s)
const Real alpha(const unsigned i) const
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
const Real getDualVariable(const unsigned i) const
unsigned maxSize_
void resetDualVariables(void)
unsigned size(void) const
const Real computeAlpha(const Real dm, const Real le) const
const Vector< Real > & subgradient(const unsigned i) const
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
unsigned remSize_
std::vector< Real > dualVariables_
Bundle_U(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)
const Real dotGi(const unsigned i, const Vector< Real > &x) const
const Real GiGj(const unsigned i, const unsigned j) const
std::vector< Real > linearizationErrors_
std::vector< Real > distanceMeasures_
virtual void initialize(const Vector< Real > &g)
void addGi(const unsigned i, const Real a, Vector< Real > &x) const
void aggregate(Vector< Real > &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const
void add(const Vector< Real > &g, const Real le, const Real dm)
const Real distanceMeasure(const unsigned i) const
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
std::vector< Ptr< Vector< Real > > > subgradients_
const Real linearizationError(const unsigned i) const
void reset(const Vector< Real > &g, const Real le, const Real dm)
void remove(const std::vector< unsigned > &ind)
void setDualVariable(const unsigned i, const Real val)
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
virtual Real dot(const Vector &x) const =0
Compute where .