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