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