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