ROL
ROL_Bundle_TT.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_TT_H
11 #define ROL_BUNDLE_TT_H
12 
13 #include "ROL_Types.hpp"
14 #include "ROL_Vector.hpp"
15 #include "ROL_StdVector.hpp"
16 #include "ROL_Bundle.hpp"
17 
18 #include "ROL_Ptr.hpp"
19 #include "ROL_LAPACK.hpp"
20 
21 #include <vector>
22 #include <limits.h>
23 #include <stdint.h>
24 #include <float.h>
25 #include <math.h>
26 #include <algorithm> // TT: std::find
27 
28 #include "ROL_LinearAlgebra.hpp"
29 
30 #define EXACT 1
31 #define TABOO_LIST 1
32 #define FIRST_VIOLATED 0
33 
40 namespace ROL {
41 
42 template<class Real>
43 class Bundle_TT : public Bundle<Real> {
44 private:
45  ROL::LAPACK<int, Real> lapack_; // TT
46 
47  int QPStatus_; // QP solver status
48  int maxind_; // maximum integer value
49  int entering_; // index of entering item
50  int LiMax_; // index of max element of diag(L)
51  int LiMin_; // index of min element of diag(L)
52 
53  unsigned maxSize_; // maximum bundle size
54  unsigned dependent_; // number of lin. dependent items in base
55  unsigned currSize_; // current size of base
56 
58  bool optimal_; // flag for optimality of restricted solution
59 
60  Real rho_;
61  Real lhNorm;
62  Real ljNorm;
63  Real lhz1_;
64  Real lhz2_;
65  Real ljz1_;
66  Real kappa_; // condition number of matrix L ( >= max|L_ii|/min|L_ii| )
67  Real objval_; // value of objective
68  Real minobjval_; // min value of objective (ever reached)
69  Real deltaLh_; // needed in case dependent row becomes independent
70  Real deltaLj_; // needed in case dependent row becomes independent
71 
72  std::vector<int> taboo_; // list of "taboo" items
73  std::vector<int> base_; // base
74 
75  LA::Matrix<Real> L_;
76  LA::Matrix<Real> Id_;
77  LA::Vector<Real> tempv_;
78  LA::Vector<Real> tempw1_;
79  LA::Vector<Real> tempw2_;
80  LA::Vector<Real> lh_;
81  LA::Vector<Real> lj_;
82  LA::Vector<Real> z1_;
83  LA::Vector<Real> z2_;
84 
85  Real GiGj(const int i, const int j) const {
87  }
88 
89  Real sgn(const Real x) const {
90  const Real zero(0), one(1);
91  return ((x < zero) ? -one :
92  ((x > zero) ? one : zero));
93  }
94 
95 public:
96  Bundle_TT(const unsigned maxSize = 10,
97  const Real coeff = 0.0,
98  const Real omega = 2.0,
99  const unsigned remSize = 2)
100  : Bundle<Real>(maxSize,coeff,omega,remSize),
101  maxSize_(maxSize), isInitialized_(false) {
102  maxind_ = std::numeric_limits<int>::max();
103  Id_.reshape(maxSize_,maxSize_);
104  for(unsigned i=0; i<maxSize_; ++i) {
105  Id_(i,i) = static_cast<Real>(1);
106  }
107  }
108 
109  unsigned solveDual(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
110  unsigned iter = 0;
111  if (Bundle<Real>::size() == 1) {
112  iter = Bundle<Real>::solveDual_dim1(t,maxit,tol);
113  }
114  else if (Bundle<Real>::size() == 2) {
115  iter = Bundle<Real>::solveDual_dim2(t,maxit,tol);
116  }
117  else {
118  iter = solveDual_arbitrary(t,maxit,tol);
119  }
120  return iter;
121  }
122 
123 /***********************************************************************************************/
124 /****************** DUAL CUTTING PLANE SUBPROBLEM ROUTINES *************************************/
125 /***********************************************************************************************/
126 private:
127  void swapRowsL(unsigned ind1, unsigned ind2, bool trans=false) {
128  const Real zero(0), one(1);
129  if( ind1 > ind2){
130  unsigned tmp = ind1;
131  ind2 = ind1;
132  ind1 = tmp;
133  }
134  unsigned dd = ind1;
135  for (unsigned n=ind1+1; n<=ind2; ++n){
136  LA::Matrix<Real> Id_n(LA::Copy,Id_,currSize_,currSize_);
137  Id_n(dd,dd) = zero; Id_n(dd,n) = one;
138  Id_n(n,dd) = one; Id_n(n,n) = zero;
139  LA::Matrix<Real> prod(currSize_,currSize_);
140  if( !trans ) {
141  prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,Id_n,L_,zero);
142  }
143  else {
144  prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,L_,Id_n,zero);
145  }
146  L_ = prod;
147  dd++;
148  }
149  }
150 
151  void updateK(void) {
152  if (currSize_ <= dependent_) { // L is empty
153  kappa_ = static_cast<Real>(1);
154  }
155  else{
156  Real tmpdiagMax = ROL_NINF<Real>();
157  Real tmpdiagMin = ROL_INF<Real>();
158  for (unsigned j=0;j<currSize_-dependent_;j++){
159  if( L_(j,j) > tmpdiagMax ){
160  tmpdiagMax = L_(j,j);
161  LiMax_ = j;
162  }
163  if( L_(j,j) < tmpdiagMin ){
164  tmpdiagMin = L_(j,j);
165  LiMin_ = j;
166  }
167  }
168  kappa_ = tmpdiagMax/tmpdiagMin;
169  }
170  }
171 
172  void addSubgradToBase(unsigned ind, Real delta) {
173  // update z1, z2, kappa
174  // swap rows if: dependent == 1 and we insert independent row (dependent row is always last)
175  // dependent == 2 and Lj has become independent (Lh still dependent)
176  if(dependent_ && (ind == currSize_-1)){
178  int tmp;
179  tmp = base_[currSize_-2];
180  base_[currSize_-2] = base_[currSize_-1];
181  base_[currSize_-1] = tmp;
182  ind--;
183  } // end if dependent
184 
185  L_(ind,ind) = delta;
186 
187  // update z1 and z2
188  unsigned zsize = ind+1;
189  z1_.resize(zsize); z2_.resize(zsize);
190  z1_[ind] = ( static_cast<Real>(1) - lhz1_ ) / delta;
191  z2_[ind] = ( Bundle<Real>::alpha(base_[ind]) - lhz2_ ) / delta;
192  //z2[zsize-1] = ( Bundle<Real>::alpha(entering_) - lhz2_ ) / delta;
193 
194  // update kappa
195  if(delta > L_(LiMax_,LiMax_)){
196  LiMax_ = ind;
197  kappa_ = delta/L_(LiMin_,LiMin_);
198  }
199  if(delta < L_(LiMin_,LiMin_)){
200  LiMin_ = ind;
201  kappa_ = L_(LiMax_,LiMax_)/delta;
202  }
203  }
204 
205  void deleteSubgradFromBase(unsigned ind, Real tol){
206  const Real zero(0), one(1);
207  // update L, currSize, base_, z1, z2, dependent, dualVariables_, kappa
208  if (ind >= currSize_-dependent_){
209  // if dependent > 0, the last one or two rows of L are lin. dependent
210  if (ind < currSize_-1){ // eliminate currSize_-2 but keep currSize_-1
211  // swap the last row with the second to last
212  swapRowsL(ind,currSize_-1);
213  base_[ind] = base_[currSize_-1];
214 #if( ! EXACT )
215  lhNorm = ljNorm; // new last row is lh
216 #endif
217  }
218 
219  dependent_--;
220  currSize_--;
221  L_.reshape(currSize_,currSize_); // the row to be eliminated is the last row
222  base_.resize(currSize_);
223 
224  // note: z1, z2, kappa need not be updated
225  return;
226  } // end if dependent item
227 
228  /* currently L_B is lower trapezoidal
229 
230  | L_1 0 0 |
231  L_B = | l d 0 |
232  | Z v L_2 |
233 
234  Apply Givens rotations to transform it to
235 
236  | L_1 0 0 |
237  | l d 0 |
238  | Z 0 L_2' |
239 
240  then delete row and column to obtain factorization of L_B' with B' = B/{i}
241 
242  L_B' = | L_1 0 |
243  | Z L_2' |
244 
245  */
246  for (unsigned j=ind+1; j<currSize_-dependent_; ++j){
247  Real ai = L_(j,ind);
248  if (std::abs(ai) <= tol*currSize_) { // nothing to do
249  continue;
250  }
251  Real aj = L_(j,j);
252  Real d, Gc, Gs;
253  // Find Givens
254  // Anderson's implementation
255  if (std::abs(aj) <= tol*currSize_){ // aj is zero
256  Gc = zero;
257  d = std::abs(ai);
258  Gs = -sgn(ai); // Gs = -sgn(ai)
259  }
260  else if ( std::abs(ai) > std::abs(aj) ){
261  Real t = aj/ai;
262  Real sgnb = sgn(ai);
263  Real u = sgnb * std::sqrt(one + t*t );
264  Gs = -one/u;
265  Gc = -Gs*t;
266  d = ai*u;
267  }
268  else{
269  Real t = ai/aj;
270  Real sgna = sgn(aj);
271  Real u = sgna * std::sqrt(one + t*t );
272  Gc = one/u;
273  Gs = -Gc*t;
274  d = aj*u;
275  }
276 
277  // // "Naive" implementation
278  // d = hypot(ai,aj);
279  // Gc = aj/d;
280  // Gs = -ai/d;
281 
282  L_(j,j) = d; L_(j,ind) = zero;
283  // apply Givens to columns i,j of L
284  for (unsigned h=j+1; h<currSize_; ++h){
285  Real tmp1 = L_(h,ind);
286  Real tmp2 = L_(h,j);
287  L_(h,ind) = Gc*tmp1 + Gs*tmp2;
288  L_(h,j) = Gc*tmp2 - Gs*tmp1;
289  }
290  // apply Givens to z1, z2
291  Real tmp1 = z1_[ind];
292  Real tmp2 = z1_[j];
293  Real tmp3 = z2_[ind];
294  Real tmp4 = z2_[j];
295  z1_[ind] = Gc*tmp1 + Gs*tmp2;
296  z1_[j] = Gc*tmp2 - Gs*tmp1;
297  z2_[ind] = Gc*tmp3 + Gs*tmp4;
298  z2_[j] = Gc*tmp4 - Gs*tmp3;
299  }// end loop over j
300 
301  if(dependent_){
302  deltaLh_ = L_(currSize_-dependent_,ind); // h = currSize_ - dependent
303  if( dependent_ > 1 ) // j = currSize_ - 1, h = currSize_ - 2
304  deltaLj_ = L_(currSize_-1,ind);
305  }
306 
307  // shift rows and columns of L by exchanging i-th row with next row and i-th column with next column until the row to be deleted is the last, then deleting last row and column
308  swapRowsL(ind,currSize_-1);
309  swapRowsL(ind,currSize_-1,true);
310  L_.reshape(currSize_-1,currSize_-1);
311 
312  // delete i-th item from z1 and z2
313  // note: z1 and z2 are of size currSize_-dependent
314  unsigned zsize = currSize_-dependent_;
315  for( unsigned m=ind; m<zsize; m++ ){
316  z1_[m] = z1_[m+1];
317  z2_[m] = z2_[m+1];
318  }
319  z1_.resize(zsize-1);
320  z2_.resize(zsize-1);
321 
322  // update base
323  base_.erase(base_.begin()+ind);
324  currSize_--; // update size
325 
326  // update kappa
327  updateK();
328 
329  if(dependent_){
330  // if some previously dependent item have become independent
331  // recompute deltaLh
332  Real ghNorm = GiGj(base_[currSize_-dependent_],base_[currSize_-dependent_]);
333  Real lhnrm(0); // exact lhNorm
334 #if( EXACT)
335  for (unsigned ii=0; ii<currSize_-dependent_; ++ii){
336  lhnrm += L_(currSize_-dependent_,ii)*L_(currSize_-dependent_,ii);
337  }
338  deltaLh_ = std::abs(ghNorm - lhnrm);
339 #else
340  Real sgn1 = sgn(deltaLh_);
341  Real sgn2 = sgn(dletaLj);
342  Real signum = sgn1 * sgn2; // sgn( deltaLh ) * sgn ( deltaLj );
343  deltaLh_ = std::abs( ghNorm - lhNorm + deltaLh_ * deltaLh_);
344 #endif
345 
346  if( std::sqrt(deltaLh_) > tol*kappa_*std::max(static_cast<Real>(1),ghNorm) ){ // originally had deltaLh without sqrt
347  unsigned newind = currSize_-dependent_;
348  dependent_--;
349  // get the last row of L
350  lh_.size(newind); // initialize to zeros;
351  lhz1_ = zero;
352  lhz2_ = zero;
353  for (unsigned ii=0; ii<newind; ++ii){
354  lh_[ii] = L_(newind,ii);
355  lhz1_ += lh_[ii]*z1_[ii];
356  lhz2_ += lh_[ii]*z2_[ii];
357  }
358  deltaLh_ = std::sqrt(deltaLh_);
359  addSubgradToBase(newind,deltaLh_);
360 
361  if(dependent_){ // dependent was 2
362 #if( ! EXACT )
363  Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
364  ljNorm -= deltaLj_ * deltaLj_;
365  lhNorm = gjNorm;
366  deltaLj_ = std::abs(gjNorm - ljNorm);
367  if ( signum < 0 )
368  deltaLj_ = - std::sqrt( deltaLj_ );
369  else
370  deltaLj_ = std::sqrt( deltaLj_ );
371 #else
372  // recompute deltaLj
373  Real gjTgh = GiGj(base_[currSize_-1],base_[currSize_-2]);
374  Real ljTlh = 0.0;
375  for (unsigned ii=0;ii<currSize_;ii++){
376  ljTlh += L_(currSize_-1,ii)*L_(currSize_-2,ii);
377  }
378  deltaLj_ = (gjTgh - ljTlh) / deltaLh_;
379 #endif
380  L_(currSize_-1,currSize_-2) = deltaLj_;
381  }
382  } // end if deltaLh > 0
383 
384  if (dependent_ > 1){ // deltaLh is 0 but deltaLj is not
385  // recompute deltaLj
386  Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
387  Real ljnrm = zero; // exact ljNorm
388 #if( EXACT )
389  for (unsigned ii=0; ii<currSize_; ++ii) {
390  ljnrm += L_(currSize_-1,ii)*L_(currSize_-1,ii);
391  }
392  deltaLj_ = std::abs(gjNorm - ljnrm);
393 #else
394  deltaLj_ = std::abs( gjNorm - ljNorm + deltaLj_ * deltaLj_);
395 #endif
396 
397  if( std::sqrt(deltaLj_) > tol*kappa_*std::max(static_cast<Real>(1),gjNorm) ){ // originally had deltaLj without sqrt
398  unsigned newind = currSize_-1;
399  dependent_--;
400  // get the last row of L
401  lj_.size(newind-1); // initialize to zeros;
402  Real ljz1__ = zero;
403  Real ljTz2 = zero;
404  for (unsigned ii=0;ii<newind-1;ii++){
405  lj_[ii] = L_(newind,ii);
406  ljz1__ += lj_[ii]*z1_[ii];
407  ljTz2 += lj_[ii]*z2_[ii];
408  }
409  deltaLj_ = std::sqrt(deltaLj_);
410  addSubgradToBase(newind,deltaLj_);
411 #if( EXACT )
412  deltaLh_ = GiGj(base_[currSize_-2],base_[currSize_-1]);
413  for (unsigned ii=0;ii<currSize_-1;ii++){
414  deltaLh_ -= L_(currSize_-2,ii)*L_(currSize_-1,ii);
415  }
416  deltaLh_ /= deltaLj_;
417 #else
418  if ( signum < 0) {
419  deltaLh_ = - std::sqrt( deltaLh_ );
420  }
421  else {
422  deltaLh_ = std::sqrt ( deltaLh_ );
423  }
424 #endif
425  L_(currSize_-1,currSize_-2) = deltaLh_;
426  } // end if deltaLj > 0
427  } // end if ( dependent > 1 )
428  } // end if(dependent)
429  }// end deleteSubgradFromBase()
430 
431  // TT: solving triangular system for TT algorithm
432  void solveSystem(int size, char tran, LA::Matrix<Real> &L, LA::Vector<Real> &v){
433  int info;
434  if( L.numRows()!=size )
435  std::cout << "Error: Wrong size matrix!" << std::endl;
436  else if( v.numRows()!=size )
437  std::cout << "Error: Wrong size vector!" << std::endl;
438  else if( size==0 )
439  return;
440  else{
441  //std::cout << L_.stride() << ' ' << size << std::endl;
442  lapack_.TRTRS( 'L', tran, 'N', size, 1, L.values(), L.stride(), v.values(), v.stride(), &info );
443  }
444  }
445 
446  // TT: check that inequality constraints are satisfied for dual variables
447  bool isFeasible(LA::Vector<Real> &v, const Real &tol){
448  bool feas = true;
449  for(int i=0;i<v.numRows();i++){
450  if(v[i]<-tol){
451  feas = false;
452  }
453  }
454  return feas;
455  }
456 
457  unsigned solveDual_TT(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
458  const Real zero(0), half(0.5), one(1);
459  Real z1z2(0), z1z1(0);
460  QPStatus_ = 1; // normal status
461  entering_ = maxind_;
462 
463  // cold start
464  optimal_ = true;
465  dependent_ = 0;
466  rho_ = ROL_INF<Real>(); // value of rho = -v
467  currSize_ = 1; // current base size
468  base_.clear();
469  base_.push_back(0); // initial base
470  L_.reshape(1,1);
471  L_(0,0) = std::sqrt(GiGj(0,0));
474  tempv_.resize(1);
475  tempw1_.resize(1);
476  tempw2_.resize(1);
477  lh_.resize(1);
478  lj_.resize(1);
479  z1_.resize(1); z2_.resize(1);
480  z1_[0] = one/L_(0,0);
481  z2_[0] = Bundle<Real>::alpha(0)/L_(0,0);
482  LiMax_ = 0; // index of max element of diag(L)
483  LiMin_ = 0; // index of min element of diag(L)
484  kappa_ = one; // condition number of matrix L ( >= max|L_ii|/min|L_ii| )
485  objval_ = ROL_INF<Real>(); // value of objective
486  minobjval_ = ROL_INF<Real>(); // min value of objective (ever reached)
487 
488  unsigned iter;
489  //-------------------------- MAIN LOOP --------------------------------//
490  for (iter=0;iter<maxit;iter++){
491  //---------------------- INNER LOOP -----------------------//
492  while( !optimal_ ){
493  switch( dependent_ ){
494  case(0): // KT system admits unique solution
495  {
496  /*
497  L = L_B'
498  */
499  z1z2 = z1_.dot(z2_);
500  z1z1 = z1_.dot(z1_);
501  rho_ = ( one + z1z2/t )/z1z1;
502  tempv_ = z1_; tempv_.scale(rho_);
503  tempw1_ = z2_; tempw1_.scale(one/t);
504  tempv_ -= tempw1_;
505  solveSystem(currSize_,'T',L_,tempv_); // tempv contains solution
506  optimal_ = true;
507  break;
508  }
509  case(1):
510  {
511  /*
512  L = | L_B' 0 | \ currSize
513  | l_h^T 0 | /
514  */
515  LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-1,currSize_-1);
516  lh_.size(currSize_-1); // initialize to zeros;
517  lhz1_ = zero;
518  lhz2_ = zero;
519  for(unsigned i=0; i<currSize_-1; ++i){
520  Real tmp = L_(currSize_-1,i);
521  lhz1_ += tmp*z1_(i);
522  lhz2_ += tmp*z2_(i);
523  lh_[i] = tmp;
524  }
525  // Test for singularity
526  if (std::abs(lhz1_-one) <= tol*kappa_){
527  // tempv is an infinite direction
528  tempv_ = lh_;
529  solveSystem(currSize_-1,'T',LBprime,tempv_);
530  tempv_.resize(currSize_); // add last entry
531  tempv_[currSize_-1] = -one;
532  optimal_ = false;
533  }
534  else{
535  // system has (unique) solution
536  rho_ = ( (Bundle<Real>::alpha(base_[currSize_-1]) - lhz2_)/t ) / ( one - lhz1_ );
537  z1z2 = z1_.dot(z2_);
538  z1z1 = z1_.dot(z1_);
539  Real tmp = ( one + z1z2 / t - rho_ * z1z1 )/( one - lhz1_ );
540  tempw1_ = z1_; tempw1_.scale(rho_);
541  tempw2_ = z2_; tempw2_.scale(one/t);
542  tempw1_ -= tempw2_;
543  tempw2_ = lh_; tempw2_.scale(tmp);
544  tempw1_ -= tempw2_;
545  solveSystem(currSize_-1,'T',LBprime,tempw1_);
546  tempv_ = tempw1_;
547  tempv_.resize(currSize_);
548  tempv_[currSize_-1] = tmp;
549  optimal_ = true;
550  }
551  break;
552  } // case(1)
553  case(2):
554  {
555  /* | L_B' 0 0 | \
556  L = | l_h^T 0 0 | | currSize
557  | l_j^T 0 0 | /
558  */
559  LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-2,currSize_-2 );
560  lj_.size(currSize_-2); // initialize to zeros;
561  lh_.size(currSize_-2); // initialize to zeros;
562  ljz1_ = zero;
563  lhz1_ = zero;
564  for(unsigned i=0; i<currSize_-2; ++i){
565  Real tmp1 = L_(currSize_-1,i);
566  Real tmp2 = L_(currSize_-2,i);
567  ljz1_ += tmp1*z1_(i);
568  lhz1_ += tmp2*z1_(i);
569  lj_[i] = tmp1;
570  lh_[i] = tmp2;
571  }
572  if(std::abs(ljz1_-one) <= tol*kappa_){
573  // tempv is an infinite direction
574  tempv_ = lj_;
575  solveSystem(currSize_-2,'T',LBprime,tempv_);
576  tempv_.resize(currSize_); // add two last entries
577  tempv_[currSize_-2] = zero;
578  tempv_[currSize_-1] = -one;
579  }
580  else{
581  // tempv is an infinite direction
582  Real mu = ( one - lhz1_ )/( one - ljz1_ );
583  tempw1_ = lj_; tempw1_.scale(-mu);
584  tempw1_ += lh_;
585  solveSystem(currSize_-2,'T',LBprime,tempw1_);
586  tempv_ = tempw1_;
587  tempv_.resize(currSize_);
588  tempv_[currSize_-2] = -one;
589  tempv_[currSize_-1] = mu;
590  }
591  optimal_ = false;
592  }// case(2)
593  } // end switch(dependent_)
594 
595  // optimal is true if tempv is a solution, otherwise, tempv is an infinite direction
596  if (optimal_){
597  // check feasibility (inequality constraints)
598  if (isFeasible(tempv_,tol*currSize_)){
599  // set dual variables to values in tempv
601  for (unsigned i=0; i<currSize_; ++i){
603  }
604  }
605  else{
606  // w_B = /bar{x}_B - x_B
607  for (unsigned i=0; i<currSize_; ++i){
609  }
610  optimal_ = false;
611  }
612  } // if(optimal)
613  else{ // choose the direction of descent
615  if ( tempv_[currSize_-1] < zero ) // w_h < 0
616  tempv_.scale(-one);
617  }
618  else{ // no i such that dualVariables_[i] == 0
619  Real sp(0);
620  for( unsigned kk=0; kk<currSize_; ++kk)
621  sp += tempv_[kk]*Bundle<Real>::alpha(base_[kk]);
622  if ( sp > zero )
623  tempv_.scale(-one);
624  }
625  }// end else ( not optimal )
626 
627  if(!optimal_){
628  // take a step in direction tempv (possibly infinite)
629  Real myeps = tol*currSize_;
630  Real step = ROL_INF<Real>();
631  for (unsigned h=0; h<currSize_; ++h){
632  if ( (tempv_[h] < -myeps) && (-Bundle<Real>::getDualVariable(base_[h])/tempv_[h] < step) )
633  if ( (Bundle<Real>::getDualVariable(base_[h]) > myeps)
634  || (Bundle<Real>::getDualVariable(base_[h]) < -myeps) ) // otherwise, consider it 0
636  }
637 
638  if (step <= zero || step == ROL_INF<Real>()){
639  QPStatus_ = -1; // invalid step
640  return iter;
641  }
642 
643  for (unsigned i=0; i<currSize_; ++i)
645  }// if(!optimal)
646 
647  //------------------------- ITEMS ELIMINATION ---------------------------//
648 
649  // Eliminate items with 0 multipliers from base
650  bool deleted = optimal_;
651  for (unsigned baseitem=0; baseitem<currSize_; ++baseitem){
652  if (Bundle<Real>::getDualVariable(base_[baseitem]) <= tol){
653  deleted = true;
654 
655 #if( TABOO_LIST )
656  // item that just entered shouldn't exit; if it does, mark it as taboo
657  if( base_[baseitem] == entering_ ){
658  taboo_.push_back(entering_);
659  entering_ = maxind_;
660  }
661 #endif
662 
663  // eliminate item from base;
664  deleteSubgradFromBase(baseitem,tol);
665 
666  } // end if(dualVariables_[baseitem] < tol)
667  } // end loop over baseitem
668 
669  if(!deleted){ // nothing deleted and not optimal
670  QPStatus_ = -2; // loop
671  return iter;
672  }
673  } // end inner loop
674 
675  Real newobjval(0), Lin(0), Quad(0); // new objective value
676  for (unsigned i=0; i<currSize_; ++i){
678  }
679  if (rho_ == ROL_NINF<Real>()){
680  Quad = -Lin/t;
681  newobjval = -half*Quad;
682  }
683  else{
684  Quad = rho_ - Lin/t;
685  newobjval = half*(rho_ + Lin/t);
686  }
687 
688 #if( TABOO_LIST )
689  // -- test for strict decrease -- //
690  // if item didn't provide decrease, move it to taboo list ...
691  if( ( entering_ < maxind_ ) && ( objval_ < ROL_INF<Real>() ) ){
692  if( newobjval >= objval_ - std::max( tol*std::abs(objval_), ROL_EPSILON<Real>() ) ){
693  taboo_.push_back(entering_);
694  }
695  }
696 #endif
697 
698  objval_ = newobjval;
699 
700  // if sufficient decrease obtained
701  if ( objval_ + std::max( tol*std::abs(objval_), ROL_EPSILON<Real>() ) <= minobjval_ ){
702  taboo_.clear(); // reset taboo list
703  minobjval_ = objval_;
704  }
705 
706  //---------------------- OPTIMALITY TEST -------------------------//
707  if ( (rho_ >= ROL_NINF<Real>()) && (objval_ <= ROL_NINF<Real>()) ) // if current x (dualVariables_) is feasible
708  break;
709 
710  entering_ = maxind_;
711  Real minro = - std::max( tol*currSize_*std::abs(objval_), ROL_EPSILON<Real>() );
712 #if ( ! FIRST_VIOLATED )
713  Real diff = ROL_NINF<Real>(), olddiff = ROL_NINF<Real>();
714 #endif
715 
716  for (unsigned bundleitem=0; bundleitem<Bundle<Real>::size(); ++bundleitem){ // loop over items in bundle
717  //for (int bundleitem=size_-1;bundleitem>-1;bundleitem--){ // loop over items in bundle (in reverse order)
718  if ( std::find(taboo_.begin(),taboo_.end(),bundleitem) != taboo_.end() ){
719  continue; // if item is taboo, move on
720  }
721 
722  if ( std::find(base_.begin(),base_.end(),bundleitem) == base_.end() ){
723  // base does not contain bundleitem
724  Real ro = zero;
725  for (unsigned j=0;j<currSize_;j++){
726  ro += Bundle<Real>::getDualVariable(base_[j]) * GiGj(bundleitem,base_[j]);
727  }
728  ro += Bundle<Real>::alpha(bundleitem)/t;
729 
730 
731  if (rho_ >= ROL_NINF<Real>()){
732  ro = ro - rho_; // note: rho = -v
733  }
734  else{
735  ro = ROL_NINF<Real>();
736  minobjval_ = ROL_INF<Real>();
737  objval_ = ROL_INF<Real>();
738  }
739 
740  if (ro < minro){
741 #if ( FIRST_VIOLATED )
742  entering_ = bundleitem;
743  break; // skip going through rest of constraints; alternatively, could look for "most violated"
744 #else
745  diff = minro - ro;
746  if ( diff > olddiff ){
747  entering_ = bundleitem;
748  olddiff = diff;
749  }
750 #endif
751  }
752 
753  } // end if item not in base
754  }// end of loop over items in bundle
755 
756  //----------------- INSERTING ITEM ------------------------//
757  if (entering_ < maxind_){ // dual constraint is violated
758  optimal_ = false;
760  base_.push_back(entering_);
761  // construct new row of L_B
762  unsigned zsize = currSize_ - dependent_; // zsize is the size of L_Bprime (current one)
763  lh_.size(zsize); // initialize to zeros;
764  lhz1_ = zero;
765  lhz2_ = zero;
766  for (unsigned i=0; i<zsize; ++i){
767  lh_[i] = GiGj(entering_,base_[i]);
768  }
769  LA::Matrix<Real> LBprime( LA::Copy,L_,zsize,zsize);
770  solveSystem(zsize,'N',LBprime,lh_); // lh = (L_B^{-1})*(G_B^T*g_h)
771  for (unsigned i=0; i<zsize; ++i){
772  lhz1_ += lh_[i]*z1_[i];
773  lhz2_ += lh_[i]*z2_[i];
774  }
775 
776  Real nrm = lh_.dot(lh_);
777  Real delta = GiGj(entering_,entering_) - nrm; // delta squared
778 #if( ! EXACT )
779  if(dependent_)
780  ljNorm = nrm; // adding second dependent
781  else
782  lhNorm = nrm; // adding first dependent
783 #endif
784 
785  currSize_++; // update base size
786 
787  L_.reshape(currSize_,currSize_);
788  zsize = currSize_ - dependent_; // zsize is the size of L_Bprime (new one)
789  for (unsigned i=0; i<zsize-1; ++i){
790  L_(currSize_-1,i) = lh_[i];
791  }
792 
793  Real deltaeps = tol*kappa_*std::max(one,lh_.dot(lh_));
794  if ( delta > deltaeps ){ // new row is independent
795  // add subgradient to the base
796  unsigned ind = currSize_-1;
797  delta = std::sqrt(delta);
798  addSubgradToBase(ind,delta);
799  }
800  else if(delta < -deltaeps){
801  dependent_++;
802  QPStatus_ = 0; // negative delta
803  return iter;
804  }
805  else{ // delta zero
806  dependent_++;
807  }
808  } // end if(entering_ < maxind_)
809  else{ // no dual constraint violated
810  if( objval_ - std::max( tol*std::abs( objval_ ), ROL_EPSILON<Real>() ) > minobjval_ ){ // check if f is as good as minf
811  QPStatus_ = -3; // loop
812  return iter;
813  }
814  }
815 
816  if(optimal_)
817  break;
818  } // end main loop
819 
820  taboo_.clear();
821  return iter;
822  }// end solveDual_TT()
823 
824  unsigned solveDual_arbitrary(const Real t, const unsigned maxit = 1000, const Real tol = 1.e-8) {
825  Real mytol = tol;
826  unsigned outermaxit = 20;
827  bool increase = false, decrease = false;
828  unsigned iter = 0;
829  for ( unsigned it=0; it < outermaxit; ++it ){
830  iter += solveDual_TT(t,maxit,mytol);
831  if ( QPStatus_ == 1 ) {
832  break;
833  }
834  else if ( QPStatus_ == -2 || QPStatus_ == -3 ) {
835  mytol /= static_cast<Real>(10);
836  decrease = true;
837  }
838  else {
839  mytol *= static_cast<Real>(10);
840  increase = true;
841  }
842  if ( (mytol > static_cast<Real>(1e-4))
843  || (mytol < static_cast<Real>(1e-16)) ){
844  break;
845  }
846  if ( increase && decrease ) {
847  break;
848  }
849  }// end outer loop
850  return iter;
851  }
852 
853 }; // class Bundle_TT
854 
855 } // namespace ROL
856 
857 #endif
858 
Bundle_TT(const unsigned maxSize=10, const Real coeff=0.0, const Real omega=2.0, const unsigned remSize=2)
unsigned currSize_
LA::Vector< Real > lh_
unsigned solveDual_TT(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
void resetDualVariables(void)
Definition: ROL_Bundle.hpp:154
bool isFeasible(LA::Vector< Real > &v, const Real &tol)
unsigned solveDual(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Real sgn(const Real x) const
unsigned solveDual_arbitrary(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
LA::Vector< Real > tempw2_
void updateK(void)
void setDualVariable(const unsigned i, const Real val)
Definition: ROL_Bundle.hpp:150
Contains definitions of custom data types in ROL.
ROL::LAPACK< int, Real > lapack_
LA::Matrix< Real > Id_
LA::Vector< Real > z1_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
LA::Vector< Real > tempw1_
const Real alpha(const unsigned i) const
Definition: ROL_Bundle.hpp:167
LA::Vector< Real > tempv_
void solveSystem(int size, char tran, LA::Matrix< Real > &L, LA::Vector< Real > &v)
std::vector< int > base_
const Real getDualVariable(const unsigned i) const
Definition: ROL_Bundle.hpp:146
LA::Vector< Real > lj_
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
unsigned dependent_
void deleteSubgradFromBase(unsigned ind, Real tol)
Provides the interface for and implements a bundle. The semidefinite quadratic subproblem is solved u...
LA::Matrix< Real > L_
void addSubgradToBase(unsigned ind, Real delta)
Real GiGj(const int i, const int j) const
std::vector< int > taboo_
unsigned solveDual_dim2(const Real t, const unsigned maxit=1000, const Real tol=1.e-8)
Definition: ROL_Bundle.hpp:299
void swapRowsL(unsigned ind1, unsigned ind2, bool trans=false)
LA::Vector< Real > z2_
Provides the interface for and implements a bundle.
Definition: ROL_Bundle.hpp:28