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