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 "Teuchos_RCP.hpp"
51 
52 #include <vector>
53 #include <set>
54 
59 namespace ROL {
60 
61 template<class Real>
62 class Bundle {
63 /***********************************************************************************************/
64 /***************** BUNDLE STORAGE **************************************************************/
65 /***********************************************************************************************/
66 private:
67  std::vector<Teuchos::RCP<Vector<Real> > > subgradients_;
68  std::vector<Real> linearizationErrors_;
69  std::vector<Real> distanceMeasures_;
70 
71  std::vector<Real> dualVariables_;
72 
73  unsigned size_;
74 
75  unsigned maxSize_;
76  unsigned remSize_;
77  Real coeff_;
78 
80 
81  void remove(const std::vector<unsigned> &ind) {
82  Real zero(0);
83  for (unsigned j = ind.back()+1; j < size_; j++) {
84  (subgradients_[j-1])->set(*(subgradients_[j]));
88  }
89  (subgradients_[size_-1])->zero();
90  linearizationErrors_[size_-1] = ROL_OVERFLOW<Real>();
91  distanceMeasures_[size_-1] = ROL_OVERFLOW<Real>();
92  dualVariables_[size_-1] = zero;
93  for (unsigned i = ind.size()-1; i > 0; --i) {
94  for (unsigned j = ind[i-1]+1; j < size_; j++) {
95  (subgradients_[j-1])->set(*(subgradients_[j]));
99  }
100  }
101  size_ -= ind.size();
102  }
103 
104  void add(const Vector<Real> &g, const Real le, const Real dm) {
105  Real zero(0);
106  (subgradients_[size_])->set(g);
108  distanceMeasures_[size_] = dm;
109  dualVariables_[size_] = zero;
110  size_++;
111  }
112 
113 protected:
114  Teuchos::RCP<Vector<Real> > tG_;
115  Teuchos::RCP<Vector<Real> > eG_;
116  Teuchos::RCP<Vector<Real> > yG_;
117 
118 /***********************************************************************************************/
119 /***************** BUNDLE MODIFICATION AND ACCESS ROUTINES *************************************/
120 /***********************************************************************************************/
121 public:
122  virtual ~Bundle(void) {}
123  Bundle(const unsigned maxSize = 10, const Real coeff = 0.0, const unsigned remSize = 2)
124  : size_(0), maxSize_(maxSize), remSize_(remSize), coeff_(coeff), isInitialized_(false) {
125  Real zero(0);
126  remSize_ = ((remSize_ < 2) ? 2 : ((remSize_ > maxSize_-1) ? maxSize_-1 : remSize_));
127  subgradients_.clear();
128  subgradients_.assign(maxSize,Teuchos::null);
129  linearizationErrors_.clear();
130  linearizationErrors_.assign(maxSize_,ROL_OVERFLOW<Real>());
131  distanceMeasures_.clear();
132  distanceMeasures_.assign(maxSize_,ROL_OVERFLOW<Real>());
133  dualVariables_.clear();
134  dualVariables_.assign(maxSize_,zero);
135  }
136 
137  void initialize(const Vector<Real> &g) {
138  if ( !isInitialized_ ) {
139  Real zero(0), one(1);
140  for (unsigned i = 0; i < maxSize_; i++) {
141  subgradients_[i] = g.clone();
142  }
143  (subgradients_[0])->set(g);
144  linearizationErrors_[0] = zero;
145  distanceMeasures_[0] = zero;
146  dualVariables_[0] = one;
147  size_++;
148  isInitialized_ = true;
149  gx_ = g.clone();
150  ge_ = g.clone();
151  tG_ = g.clone();
152  yG_ = g.clone();
153  eG_ = g.clone();
154  }
155  }
156 
157  const Real linearizationError(const unsigned i) const {
158  return linearizationErrors_[i];
159  }
160 
161  const Real distanceMeasure(const unsigned i) const {
162  return distanceMeasures_[i];
163  }
164 
165  const Vector<Real> & subgradient(const unsigned i) const {
166  return *(subgradients_[i]);
167  }
168 
169  const Real computeAlpha(const Real dm, const Real le) const {
170  Real alpha = le, two(2);
171  if ( coeff_ > ROL_EPSILON<Real>() ) {
172  alpha = std::max(coeff_*std::pow(dm,two),le);
173  }
174  return alpha;
175  }
176 
177  const Real alpha(const unsigned i) const {
179  }
180 
181  unsigned size(void) const {
182  return size_;
183  }
184 
185  void aggregate(Vector<Real> &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const {
186  Real zero(0), one(1);
187  aggSubGrad.zero(); aggLinErr = zero; aggDistMeas = zero; eG_->zero();
188  Real eLE(0), eDM(0), yLE(0), yDM(0), tLE(0), tDM(0);
189  for (unsigned i = 0; i < size_; i++) {
190  // Compute aggregate subgradient using Kahan's compensated sum
191  tG_->set(aggSubGrad);
192  yG_->set(*subgradients_[i]); yG_->scale(dualVariables_[i]); yG_->plus(*eG_);
193  aggSubGrad.set(*tG_); aggSubGrad.plus(*yG_);
194  eG_->set(*tG_); eG_->axpy(-one,aggSubGrad); eG_->plus(*yG_);
195  // Compute aggregate linearization error using Kahan's compensated sum
196  tLE = aggLinErr;
197  yLE = dualVariables_[i]*linearizationErrors_[i] + eLE;
198  aggLinErr = tLE + yLE;
199  eLE = (tLE - aggLinErr) + yLE;
200  // Compute aggregate distance measure using Kahan's compensated sum
201  tDM = aggDistMeas;
202  yDM = dualVariables_[i]*distanceMeasures_[i] + eDM;
203  aggDistMeas = tDM + yDM;
204  eDM = (tDM - aggDistMeas) + yDM;
205  }
206  }
207 
208  void reset(const Vector<Real> &g, const Real le, const Real dm) {
209  if (size_ == maxSize_) {
210  // Find indices to remove
211  unsigned loc = size_, cnt = 0;
212  std::vector<unsigned> ind(remSize_,0);
213  for (unsigned i = size_; i > 0; --i) {
214  if ( std::abs(linearizationErrors_[i-1]) < ROL_EPSILON<Real>() ) {
215  loc = i-1;
216  break;
217  }
218  }
219  for (unsigned i = 0; i < size_; i++) {
220  if ( i < loc || i > loc ) {
221  ind[cnt] = i;
222  cnt++;
223  }
224  if (cnt == remSize_) {
225  break;
226  }
227  }
228  // Remove indices
229  remove(ind);
230  // Add aggregate subgradient
231  add(g,le,dm);
232  }
233  }
234 
235  void update(const bool flag, const Real linErr, const Real distMeas,
236  const Vector<Real> &g, const Vector<Real> &s) {
237  Real zero(0);
238  if ( flag ) {
239  // Serious step taken: Update linearlization errors and distance measures
240  for (unsigned i = 0; i < size_; i++) {
241  linearizationErrors_[i] += linErr - subgradients_[i]->dot(s.dual());
242  distanceMeasures_[i] += distMeas;
243  }
244  linearizationErrors_[size_] = zero;
245  distanceMeasures_[size_] = zero;
246  }
247  else {
248  // Null step taken:
249  linearizationErrors_[size_] = linErr;
250  distanceMeasures_[size_] = distMeas;
251  }
252  // Update (sub)gradient bundle
253  (subgradients_[size_])->set(g);
254  // Update dual variables
255  dualVariables_[size_] = zero;
256  // Update bundle size
257  size_++;
258  }
259 
260  // TT: adding access routines for derived classes
261 protected:
262 
263  Real getDualVariables (const unsigned i){
264  return dualVariables_[i];
265  }
266 
267  void setDualVariables(const unsigned i, const Real val) {
268  dualVariables_[i] = val;
269  }
270 
271  void resetDualVariables(void){
272  Real zero(0);
273  dualVariables_.assign(size_,zero);
274  }
275 
276 /***********************************************************************************************/
277 /***************** DUAL CUTTING PLANE PROBLEM ROUTINES *****************************************/
278 /***********************************************************************************************/
279 protected:
280  Teuchos::RCP<Vector<Real> > gx_;
281  Teuchos::RCP<Vector<Real> > ge_;
282 
283 private:
284  std::set<unsigned> workingSet_;
285  std::set<unsigned> nworkingSet_;
286 
287  void initializeDualSolver(void) {
288 // for (unsigned i = 0; i < maxSize_; i++) {
289 // dualVariables_[i] = 0.0;
290 // }
291 // dualVariables_[0] = 1.0;
292 // for (unsigned i = 0; i < maxSize_; i++) {
293 // dualVariables_[i] = ((i<size_) ? 1.0/(Real)size_ : 0.0);
294 // }
295 // nworkingSet_.clear();
296 // workingSet_.clear();
297 // for (unsigned i = 0; i < size_; i++) {
298 // nworkingSet_.insert(i);
299 // }
300  Real sum(0), err(0), tmp(0), y(0), zero(0);
301  for (unsigned i = 0; i < size_; i++) {
302  // Compute sum of dualVariables_ using Kahan's compensated sum
303  tmp = sum;
304  y = dualVariables_[i] + err;
305  sum = tmp + y;
306  err = (tmp - sum) + y;
307  }
308  for (unsigned i = 0; i < size_; i++) {
309  dualVariables_[i] /= sum;
310  }
311  nworkingSet_.clear();
312  workingSet_.clear();
313  for (unsigned i = 0; i < size_; i++) {
314  if ( dualVariables_[i] == zero ) {
315  workingSet_.insert(i);
316  }
317  else {
318  nworkingSet_.insert(i);
319  }
320  }
321  }
322 
323  Real evaluateObjective(std::vector<Real> &g, const std::vector<Real> &x, const Real t) const {
324  Real one(1), half(0.5);
325  gx_->zero(); eG_->zero();
326  for (unsigned i = 0; i < size_; i++) {
327  // Compute Gx using Kahan's compensated sum
328  tG_->set(*gx_);
329  yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[i]));
330  gx_->set(*tG_); gx_->plus(*yG_);
331  eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->plus(*yG_);
332  }
333  Real Hx(0), val(0), err(0), tmp(0), y(0);
334  for (unsigned i = 0; i < size_; i++) {
335  // Compute < g_i, Gx >
336  Hx = gx_->dot(*(subgradients_[i]));
337  // Add to the objective function value using Kahan's compensated sum
338  tmp = val;
339  y = x[i]*(half*Hx + alpha(i)/t) + err;
340  val = tmp + y;
341  err = (tmp - val) + y;
342  // Add gradient component
343  g[i] = Hx + alpha(i)/t;
344  }
345  return val;
346  }
347 
348  void applyFullMatrix(std::vector<Real> &Hx, const std::vector<Real> &x) const {
349  Real one(1);
350  gx_->zero(); eG_->zero();
351  for (unsigned i = 0; i < size_; i++) {
352  // Compute Gx using Kahan's compensated sum
353  tG_->set(*gx_);
354  yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[i]));
355  gx_->set(*tG_); gx_->plus(*yG_);
356  eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->plus(*yG_);
357  }
358  for (unsigned i = 0; i < size_; i++) {
359  // Compute < g_i, Gx >
360  Hx[i] = subgradients_[i]->dot(*gx_);
361  }
362  }
363 
364  void applyMatrix(std::vector<Real> &Hx, const std::vector<Real> &x) const {
365  Real one(1);
366  gx_->zero(); eG_->zero();
367  unsigned n = nworkingSet_.size();
368  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
369  for (unsigned i = 0; i < n; i++) {
370  // Compute Gx using Kahan's compensated sum
371  tG_->set(*gx_);
372  yG_->set(*eG_); yG_->axpy(x[i],*(subgradients_[*it]));
373  gx_->set(*tG_); gx_->plus(*yG_);
374  eG_->set(*tG_); eG_->axpy(-one,*gx_); eG_->plus(*yG_);
375  it++;
376  }
377  it = nworkingSet_.begin();
378  for (unsigned i = 0; i < n; i++) {
379  // Compute < g_i, Gx >
380  Hx[i] = subgradients_[*it]->dot(*gx_); it++;
381  }
382  }
383 
384  void computeLagMult(std::vector<Real> &lam, const Real mu, const std::vector<Real> &g) const {
385  Real zero(0);
386  unsigned n = workingSet_.size();
387  if ( n > 0 ) {
388  lam.resize(n,zero);
389  typename std::set<unsigned>::iterator it = workingSet_.begin();
390  for (unsigned i = 0; i < n; i++) {
391  lam[i] = g[*it] - mu; it++;
392  }
393  }
394  else {
395  lam.clear();
396  }
397  }
398 
399  bool isNonnegative(unsigned &ind, const std::vector<Real> &x) const {
400  bool flag = true;
401  unsigned n = workingSet_.size(); ind = size_;
402  if ( n > 0 ) {
403  Real min = ROL_OVERFLOW<Real>();
404  typename std::set<unsigned>::iterator it = workingSet_.begin();
405  for (unsigned i = 0; i < n; i++) {
406  if ( x[i] < min ) {
407  ind = *it;
408  min = x[i];
409  }
410  it++;
411  }
412  flag = ((min < -ROL_EPSILON<Real>()) ? false : true);
413  }
414  return flag;
415  }
416 
417  Real computeAlpha(unsigned &ind, const std::vector<Real> &x, const std::vector<Real> &p) const {
418  Real alpha(1), tmp(0), zero(0); ind = size_;
419  typename std::set<unsigned>::iterator it;
420  for (it = nworkingSet_.begin(); it != nworkingSet_.end(); it++) {
421  if ( p[*it] < -ROL_EPSILON<Real>() ) {
422  tmp = -x[*it]/p[*it];
423  if ( alpha >= tmp ) {
424  alpha = tmp;
425  ind = *it;
426  }
427  }
428  }
429  return std::max(zero,alpha);
430  }
431 
432  unsigned solveEQPsubproblem(std::vector<Real> &s, Real &mu,
433  const std::vector<Real> &g, const Real tol) const {
434  // Build reduced QP information
435  Real zero(0);
436  unsigned n = nworkingSet_.size(), cnt = 0;
437  mu = zero;
438  s.assign(size_,zero);
439  if ( n > 0 ) {
440  std::vector<Real> gk(n,zero);
441  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
442  for (unsigned i = 0; i < n; i++) {
443  gk[i] = g[*it]; it++;
444  }
445  std::vector<Real> sk(n,zero);
446  cnt = projectedCG(sk,mu,gk,tol);
447  it = nworkingSet_.begin();
448  for (unsigned i = 0; i < n; i++) {
449  s[*it] = sk[i]; it++;
450  }
451  }
452  return cnt;
453  }
454 
455  void applyPreconditioner(std::vector<Real> &Px, const std::vector<Real> &x) const {
456  Real zero(0);
457  int type = 0;
458  std::vector<Real> tmp(Px.size(),zero);
459  switch (type) {
460  case 0: applyPreconditioner_Identity(tmp,x); break;
461  case 1: applyPreconditioner_Jacobi(tmp,x); break;
462  case 2: applyPreconditioner_SymGS(tmp,x); break;
463  }
465  }
466 
467  void applyG(std::vector<Real> &Gx, const std::vector<Real> &x) const {
468  int type = 0;
469  switch (type) {
470  case 0: applyG_Identity(Gx,x); break;
471  case 1: applyG_Jacobi(Gx,x); break;
472  case 2: applyG_SymGS(Gx,x); break;
473  }
474  }
475 
476  void applyPreconditioner_Identity(std::vector<Real> &Px, const std::vector<Real> &x) const {
477  unsigned dim = nworkingSet_.size();
478  Real sum(0), err(0), tmp(0), y(0);
479  for (unsigned i = 0; i < dim; i++) {
480  // Compute sum of x using Kahan's compensated sum
481  tmp = sum;
482  y = x[i] + err;
483  sum = tmp + y;
484  err = (tmp - sum) + y;
485  }
486  sum /= (Real)dim;
487  for (unsigned i = 0; i < dim; i++) {
488  Px[i] = x[i] - sum;
489  }
490  }
491 
492  void applyG_Identity(std::vector<Real> &Gx, const std::vector<Real> &x) const {
493  Gx.assign(x.begin(),x.end());
494  }
495 
496  void applyPreconditioner_Jacobi(std::vector<Real> &Px, const std::vector<Real> &x) const {
497  unsigned dim = nworkingSet_.size();
498  Real eHe(0), sum(0), one(1), zero(0);
499  Real errX(0), tmpX(0), yX(0), errE(0), tmpE(0), yE(0);
500  std::vector<Real> gg(dim,zero);
501  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
502  for (unsigned i = 0; i < dim; i++) {
503  gg[i] = one/std::abs(subgradients_[*it]->dot(*(subgradients_[*it]))); it++;
504  // Compute sum of inv(D)x using Kahan's aggregated sum
505  tmpX = sum;
506  yX = x[i]*gg[i] + errX;
507  sum = tmpX + errX;
508  errX = (tmpX - sum) + yX;
509  // Compute sum of inv(D)e using Kahan's aggregated sum
510  tmpE = eHe;
511  yE = gg[i] + errE;
512  eHe = tmpE + yE;
513  errE = (tmpE - eHe) + yE;
514  }
515  sum /= eHe;
516  for (unsigned i = 0; i < dim; i++) {
517  Px[i] = (x[i]-sum)*gg[i];
518  }
519  }
520 
521  void applyG_Jacobi(std::vector<Real> &Gx, const std::vector<Real> &x) const {
522  unsigned dim = nworkingSet_.size();
523  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
524  for (unsigned i = 0; i < dim; i++) {
525  Gx[i] = std::abs(subgradients_[*it]->dot(*(subgradients_[*it])))*x[i]; it++;
526  }
527  }
528 
529  void applyPreconditioner_SymGS(std::vector<Real> &Px, const std::vector<Real> &x) const {
530  int dim = nworkingSet_.size();
531  //unsigned cnt = 0;
532  gx_->zero(); ge_->zero();
533  Real eHx(0), eHe(0), one(1);
534  // Forward substitution
535  std::vector<Real> x1(dim,0), e1(dim,0),gg(dim,0);
536  typename std::set<unsigned>::iterator it, jt;
537  it = nworkingSet_.begin();
538  for (int i = 0; i < dim; i++) {
539  gx_->zero(); ge_->zero(); jt = nworkingSet_.begin();
540  for (int j = 0; j < i; j++) {
541  gx_->axpy(x1[j],*(subgradients_[*jt]));
542  ge_->axpy(e1[j],*(subgradients_[*jt]));
543  jt++;
544  }
545  gg[i] = subgradients_[*it]->dot(*(subgradients_[*it]));
546  x1[i] = (x[i] - gx_->dot(*(subgradients_[*it])))/gg[i];
547  e1[i] = (one - ge_->dot(*(subgradients_[*it])))/gg[i];
548  it++;
549  }
550  // Apply diagonal
551  for (int i = 0; i < dim; i++) {
552  x1[i] *= gg[i];
553  e1[i] *= gg[i];
554  }
555  // Back substitution
556  std::vector<Real> Hx(dim,0), He(dim,0); it = nworkingSet_.end();
557  for (int i = dim-1; i >= 0; --i) {
558  it--;
559  gx_->zero(); ge_->zero(); jt = nworkingSet_.end();
560  for (int j = dim-1; j >= i+1; --j) {
561  jt--;
562  gx_->axpy(Hx[j],*(subgradients_[*jt]));
563  ge_->axpy(He[j],*(subgradients_[*jt]));
564  }
565  Hx[i] = (x1[i] - gx_->dot(*(subgradients_[*it])))/gg[i];
566  He[i] = (e1[i] - ge_->dot(*(subgradients_[*it])))/gg[i];
567  // Compute sums
568  eHx += Hx[i];
569  eHe += He[i];
570  }
571  // Accumulate preconditioned vector
572  for (int i = 0; i < dim; i++) {
573  Px[i] = Hx[i] - (eHx/eHe)*He[i];
574  }
575  }
576 
577  void applyG_SymGS(std::vector<Real> &Gx, const std::vector<Real> &x) const {
578  unsigned dim = nworkingSet_.size();
579  typename std::set<unsigned>::iterator it = nworkingSet_.begin();
580  for (unsigned i = 0; i < dim; i++) {
581  Gx[i] = std::abs(subgradients_[*it]->dot(*(subgradients_[*it])))*x[i]; it++;
582  }
583  }
584 
585  void computeResidualUpdate(std::vector<Real> &r, std::vector<Real> &g) const {
586  unsigned n = g.size();
587  std::vector<Real> Gg(n,0);
588  Real y(0), ytmp(0), yprt(0), yerr(0);
589  applyPreconditioner(g,r);
590  applyG(Gg,g);
591  // Compute multiplier using Kahan's compensated sum
592  for (unsigned i = 0; i < n; i++) {
593  ytmp = y;
594  yprt = (r[i] - Gg[i]) + yerr;
595  y = ytmp + yprt;
596  yerr = (ytmp - y) + yprt;
597  }
598  y /= (Real)n;
599  for (unsigned i = 0; i < n; i++) {
600  r[i] -= y;
601  }
602  applyPreconditioner(g,r);
603  }
604 
605  unsigned projectedCG(std::vector<Real> &x, Real &mu, const std::vector<Real> &b, const Real tol) const {
606  Real one(1), zero(0);
607  unsigned n = nworkingSet_.size();
608  std::vector<Real> r(n,0), r0(n,0), g(n,0), d(n,0), Ad(n,0);
609  // Compute residual Hx + g = g with x = 0
610  x.assign(n,0);
611  scale(r,one,b);
612  r0.assign(r.begin(),r.end());
613  // Precondition residual
615  Real rg = dot(r,g), rg0(0);
616  // Get search direction
617  scale(d,-one,g);
618  Real alpha(0), kappa(0), beta(0), TOL(1.e-2);
619  Real CGtol = std::min(tol,TOL*rg);
620  unsigned cnt = 0;
621  while (rg > CGtol && cnt < 2*n+1) {
622  applyMatrix(Ad,d);
623  kappa = dot(d,Ad);
624  alpha = rg/kappa;
625  axpy(alpha,d,x);
626  axpy(alpha,Ad,r);
627  axpy(alpha,Ad,r0);
629  rg0 = rg;
630  rg = dot(r,g);
631  beta = rg/rg0;
632  scale(d,beta);
633  axpy(-one,g,d);
634  cnt++;
635  }
636  // Compute multiplier for equality constraint using Kahan's compensated sum
637  mu = zero;
638  Real err(0), tmp(0), y(0);
639  for (unsigned i = 0; i < n; i++) {
640  tmp = mu;
641  y = r0[i] + err;
642  mu = tmp + y;
643  err = (tmp - mu) + y;
644  }
645  mu /= (Real)n;
646  // Return iteration count
647  return cnt;
648  }
649 
650  Real dot(const std::vector<Real> &x, const std::vector<Real> &y) const {
651  // Compute dot product of two vectors using Kahan's compensated sum
652  Real val(0), err(0), tmp(0), y0(0);
653  unsigned n = std::min(x.size(),y.size());
654  for (unsigned i = 0; i < n; i++) {
655  tmp = val;
656  y0 = x[i]*y[i] + err;
657  val = tmp + y0;
658  err = (tmp - val) + y0;
659  }
660  return val;
661  }
662 
663  Real norm(const std::vector<Real> &x) const {
664  return std::sqrt(dot(x,x));
665  }
666 
667  void axpy(const Real a, const std::vector<Real> &x, std::vector<Real> &y) const {
668  unsigned n = std::min(y.size(),x.size());
669  for (unsigned i = 0; i < n; i++) {
670  y[i] += a*x[i];
671  }
672  }
673 
674  void scale(std::vector<Real> &x, const Real a) const {
675  for (unsigned i = 0; i < x.size(); i++) {
676  x[i] *= a;
677  }
678  }
679 
680  void scale(std::vector<Real> &x, const Real a, const std::vector<Real> &y) const {
681  unsigned n = std::min(x.size(),y.size());
682  for (unsigned i = 0; i < n; i++) {
683  x[i] = a*y[i];
684  }
685  }
686 
687  unsigned solveDual_dim1(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
688  dualVariables_[0] = (Real)1;
689  //std::cout << "dim = " << size_ << " iter = " << 0 << " CONVERGED!\n";
690  return 0;
691  }
692 
693  unsigned solveDual_dim2(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
694  Real diffg = gx_->dot(*gx_), zero(0), one(1), half(0.5);
695  gx_->set(*subgradients_[0]); gx_->axpy(-one,*subgradients_[1]);
696  if ( std::abs(diffg) > ROL_EPSILON<Real>() ) {
697  Real diffa = (alpha(0)-alpha(1))/t;
698  Real gdiffg = subgradients_[1]->dot(*gx_);
699  dualVariables_[0] = std::min(one,std::max(zero,-(gdiffg+diffa)/diffg));
700  dualVariables_[1] = one-dualVariables_[0];
701  }
702  else {
703  if ( std::abs(alpha(0)-alpha(1)) > ROL_EPSILON<Real>() ) {
704  if ( alpha(0) < alpha(1) ) {
705  dualVariables_[0] = one; dualVariables_[1] = zero;
706  }
707  else if ( alpha(0) > alpha(1) ) {
708  dualVariables_[0] = zero; dualVariables_[1] = one;
709  }
710  }
711  else {
712  dualVariables_[0] = half; dualVariables_[1] = half;
713  }
714  }
715  //std::cout << "dim = " << size_ << " iter = " << 0 << " CONVERGED!\n";
716  return 0;
717  }
718 
719  unsigned solveDual_arbitrary(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
721  bool nonneg = false;
722  unsigned ind = 0, i = 0, CGiter = 0;
723  Real snorm(0), alpha(0), mu(0), one(1), zero(0);
724  std::vector<Real> s(size_,0), Hs(size_,0), g(size_,0), lam(size_+1,0);
725  //Real val = evaluateObjective(g,dualVariables_,t);
727  for (i = 0; i < maxit; i++) {
728  CGiter += solveEQPsubproblem(s,mu,g,tol);
729  snorm = norm(s);
730  if ( snorm < ROL_EPSILON<Real>() ) {
731  computeLagMult(lam,mu,g);
732  nonneg = isNonnegative(ind,lam);
733  if ( nonneg ) {
734  break;
735  }
736  else {
737  alpha = one;
738  if ( ind < size_ ) {
739  workingSet_.erase(ind);
740  nworkingSet_.insert(ind);
741  }
742  }
743  }
744  else {
746  if ( alpha > zero ) {
748  applyFullMatrix(Hs,s);
749  axpy(alpha,Hs,g);
750  }
751  if (ind < size_) {
752  workingSet_.insert(ind);
753  nworkingSet_.erase(ind);
754  }
755  }
756  //std::cout << "iter = " << i << " snorm = " << snorm << " alpha = " << alpha << "\n";
757  }
758  //Real crit = computeCriticality(g);
759  //std::cout << "Criticality Measure: " << crit << "\n";
760  //std::cout << "dim = " << size_ << " iter = " << i << " CGiter = " << CGiter << " CONVERGED!\n";
761  return i;
762  }
763 
764 public:
765  virtual unsigned solveDual(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
766  unsigned iter = 0;
767  if (size_ == 1) {
768  iter = solveDual_dim1(t,maxit,tol);
769  }
770  else if (size_ == 2) {
771  iter = solveDual_dim2(t,maxit,tol);
772  }
773  else {
774  iter = solveDual_arbitrary(t,maxit,tol);
775  }
776  return iter;
777  }
778 
779 private:
780  /************************************************************************/
781  /********************** PROJECTION ONTO FEASIBLE SET ********************/
782  /************************************************************************/
783  void project(std::vector<Real> &x, const std::vector<Real> &v) const {
784  std::vector<Real> vsort(size_,0);
785  vsort.assign(v.begin(),v.end());
786  std::sort(vsort.begin(),vsort.end());
787  Real sum(-1), lam(0), zero(0), one(1);
788  for (int i = size_-1; i > 0; i--) {
789  sum += vsort[i];
790  if ( sum >= ((Real)(size_-i))*vsort[i-1] ) {
791  lam = sum/(Real)(size_-i);
792  break;
793  }
794  }
795  if (lam == zero) {
796  lam = (sum+vsort[0])/(Real)size_;
797  }
798  for (int i = 0; i < size_; i++) {
799  x[i] = std::max(zero,v[i] - lam);
800  }
801  }
802 
803  Real computeCriticality(const std::vector<Real> &g) {
804  Real zero(0), one(1);
805  std::vector<Real> x(size_,0), Px(size_,0);
806  axpy(one,dualVariables_,x);
807  axpy(-one,g,x);
808  project(Px,x);
809  scale(x,zero);
810  axpy(one,dualVariables_,x);
811  axpy(-one,Px,x);
812  return norm(x);
813  }
814 }; // class Bundle
815 
816 } // namespace ROL
817 
818 #endif
819 
820 // void aggregate(Vector<Real> &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const {
821 // aggSubGrad.zero(); aggLinErr = 0.0; aggDistMeas = 0.0;
822 // for (unsigned i = 0; i < size_; i++) {
823 // //if ( dualVariables_[i] > ROL_EPSILON<Real>() ) {
824 // aggSubGrad.axpy(dualVariables_[i],*(subgradients_[i]));
825 // aggLinErr += dualVariables_[i]*linearizationErrors_[i];
826 // aggDistMeas += dualVariables_[i]*distanceMeasures_[i];
827 // //}
828 // }
829 // }
unsigned maxSize_
Definition: ROL_Bundle.hpp:75
void initializeDualSolver(void)
Definition: ROL_Bundle.hpp:287
void applyG_Jacobi(std::vector< Real > &Gx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:521
void resetDualVariables(void)
Definition: ROL_Bundle.hpp:271
void applyG(std::vector< Real > &Gx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:467
unsigned projectedCG(std::vector< Real > &x, Real &mu, const std::vector< Real > &b, const Real tol) const
Definition: ROL_Bundle.hpp:605
void computeResidualUpdate(std::vector< Real > &r, std::vector< Real > &g) const
Definition: ROL_Bundle.hpp:585
void applyFullMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:348
std::set< unsigned > workingSet_
Definition: ROL_Bundle.hpp:284
virtual void plus(const Vector &x)=0
Compute , where .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y) const
Definition: ROL_Bundle.hpp:650
void update(const bool flag, const Real linErr, const Real distMeas, const Vector< Real > &g, const Vector< Real > &s)
Definition: ROL_Bundle.hpp:235
unsigned size_
Definition: ROL_Bundle.hpp:73
Teuchos::RCP< Vector< Real > > yG_
Definition: ROL_Bundle.hpp:116
const Real alpha(const unsigned i) const
Definition: ROL_Bundle.hpp:177
Contains definitions of custom data types in ROL.
Real computeCriticality(const std::vector< Real > &g)
Definition: ROL_Bundle.hpp:803
const Vector< Real > & subgradient(const unsigned i) const
Definition: ROL_Bundle.hpp:165
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void reset(const Vector< Real > &g, const Real le, const Real dm)
Definition: ROL_Bundle.hpp:208
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void applyPreconditioner_SymGS(std::vector< Real > &Px, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:529
std::set< unsigned > nworkingSet_
Definition: ROL_Bundle.hpp:285
std::vector< Real > dualVariables_
Definition: ROL_Bundle.hpp:71
const Real distanceMeasure(const unsigned i) const
Definition: ROL_Bundle.hpp:161
Teuchos::RCP< Vector< Real > > ge_
Definition: ROL_Bundle.hpp:281
void applyG_SymGS(std::vector< Real > &Gx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:577
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:213
void scale(std::vector< Real > &x, const Real a, const std::vector< Real > &y) const
Definition: ROL_Bundle.hpp:680
std::vector< Teuchos::RCP< Vector< Real > > > subgradients_
Definition: ROL_Bundle.hpp:67
void computeLagMult(std::vector< Real > &lam, const Real mu, const std::vector< Real > &g) const
Definition: ROL_Bundle.hpp:384
virtual ~Bundle(void)
Definition: ROL_Bundle.hpp:122
Real norm(const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:663
unsigned remSize_
Definition: ROL_Bundle.hpp:76
Real computeAlpha(unsigned &ind, const std::vector< Real > &x, const std::vector< Real > &p) const
Definition: ROL_Bundle.hpp:417
Real getDualVariables(const unsigned i)
Definition: ROL_Bundle.hpp:263
std::vector< Real > distanceMeasures_
Definition: ROL_Bundle.hpp:69
void applyPreconditioner_Jacobi(std::vector< Real > &Px, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:496
unsigned size(void) const
Definition: ROL_Bundle.hpp:181
const Real computeAlpha(const Real dm, const Real le) const
Definition: ROL_Bundle.hpp:169
void applyPreconditioner(std::vector< Real > &Px, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:455
unsigned solveDual_dim1(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Definition: ROL_Bundle.hpp:687
void axpy(const Real a, const std::vector< Real > &x, std::vector< Real > &y) const
Definition: ROL_Bundle.hpp:667
Teuchos::RCP< Vector< Real > > eG_
Definition: ROL_Bundle.hpp:115
Teuchos::RCP< Vector< Real > > tG_
Definition: ROL_Bundle.hpp:114
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void applyPreconditioner_Identity(std::vector< Real > &Px, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:476
void applyG_Identity(std::vector< Real > &Gx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:492
void applyMatrix(std::vector< Real > &Hx, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:364
unsigned solveEQPsubproblem(std::vector< Real > &s, Real &mu, const std::vector< Real > &g, const Real tol) const
Definition: ROL_Bundle.hpp:432
void setDualVariables(const unsigned i, const Real val)
Definition: ROL_Bundle.hpp:267
virtual unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Definition: ROL_Bundle.hpp:765
std::vector< Real > linearizationErrors_
Definition: ROL_Bundle.hpp:68
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Definition: ROL_Bundle.hpp:693
void initialize(const Vector< Real > &g)
Definition: ROL_Bundle.hpp:137
void scale(std::vector< Real > &x, const Real a) const
Definition: ROL_Bundle.hpp:674
void add(const Vector< Real > &g, const Real le, const Real dm)
Definition: ROL_Bundle.hpp:104
Real evaluateObjective(std::vector< Real > &g, const std::vector< Real > &x, const Real t) const
Definition: ROL_Bundle.hpp:323
Teuchos::RCP< Vector< Real > > gx_
Definition: ROL_Bundle.hpp:280
bool isInitialized_
Definition: ROL_Bundle.hpp:79
Bundle(const unsigned maxSize=10, const Real coeff=0.0, const unsigned remSize=2)
Definition: ROL_Bundle.hpp:123
bool isNonnegative(unsigned &ind, const std::vector< Real > &x) const
Definition: ROL_Bundle.hpp:399
const Real linearizationError(const unsigned i) const
Definition: ROL_Bundle.hpp:157
Provides the interface for and implments a bundle.
Definition: ROL_Bundle.hpp:62
void project(std::vector< Real > &x, const std::vector< Real > &v) const
Definition: ROL_Bundle.hpp:783
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Definition: ROL_Bundle.hpp:719
void aggregate(Vector< Real > &aggSubGrad, Real &aggLinErr, Real &aggDistMeas) const
Definition: ROL_Bundle.hpp:185