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