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