10 #ifndef ROL_BUNDLE_U_TT_DEF_H
11 #define ROL_BUNDLE_U_TT_DEF_H
22 #define FIRST_VIOLATED 0
26 template<
typename Real>
30 const unsigned remSize)
31 :
Bundle_U<Real>(maxSize,coeff,omega,remSize),
32 maxSize_(maxSize), isInitialized_(false) {
33 maxind_ = std::numeric_limits<int>::max();
36 Id_(i,i) =
static_cast<Real
>(1);
40 template<
typename Real>
42 const Real
zero(0), one(1);
43 return ((x <
zero) ? -one :
47 template<
typename Real>
57 iter = solveDual_arbitrary(t,maxit,tol);
62 template<
typename Real>
64 const Real
zero(0), one(1);
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_);
77 prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,Id_n,L_,
zero);
80 prod.multiply(LA::ETransp::NO_TRANS,LA::ETransp::NO_TRANS,one,L_,Id_n,
zero);
87 template<
typename Real>
89 if (currSize_ <= dependent_) {
90 kappa_ =
static_cast<Real
>(1);
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 ){
100 if( L_(j,j) < tmpdiagMin ){
101 tmpdiagMin = L_(j,j);
105 kappa_ = tmpdiagMax/tmpdiagMin;
109 template<
typename Real>
114 if(dependent_ && (ind == currSize_-1)){
115 swapRowsL(currSize_-2,currSize_-1);
117 tmp = base_[currSize_-2];
118 base_[currSize_-2] = base_[currSize_-1];
119 base_[currSize_-1] = tmp;
126 unsigned zsize = ind+1;
127 z1_.resize(zsize); z2_.resize(zsize);
128 z1_[ind] = (
static_cast<Real
>(1) - lhz1_ ) / delta;
133 if(delta > L_(LiMax_,LiMax_)){
135 kappa_ = delta/L_(LiMin_,LiMin_);
137 if(delta < L_(LiMin_,LiMin_)){
139 kappa_ = L_(LiMax_,LiMax_)/delta;
143 template<
typename Real>
145 const Real
zero(0), one(1);
147 if (ind >= currSize_-dependent_){
149 if (ind < currSize_-1){
151 swapRowsL(ind,currSize_-1);
152 base_[ind] = base_[currSize_-1];
160 L_.reshape(currSize_,currSize_);
161 base_.resize(currSize_);
185 for (
unsigned j=ind+1; j<currSize_-dependent_; ++j){
187 if (std::abs(ai) <= tol*currSize_) {
194 if (std::abs(aj) <= tol*currSize_){
199 else if ( std::abs(ai) > std::abs(aj) ){
202 Real u = sgnb * std::sqrt(one + t*t );
210 Real u = sgna * std::sqrt(one + t*t );
221 L_(j,j) = d; L_(j,ind) =
zero;
223 for (
unsigned h=j+1; h<currSize_; ++h){
224 Real tmp1 = L_(h,ind);
226 L_(h,ind) = Gc*tmp1 + Gs*tmp2;
227 L_(h,j) = Gc*tmp2 - Gs*tmp1;
230 Real tmp1 = z1_[ind];
232 Real tmp3 = z2_[ind];
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;
241 deltaLh_ = L_(currSize_-dependent_,ind);
243 deltaLj_ = L_(currSize_-1,ind);
247 swapRowsL(ind,currSize_-1);
248 swapRowsL(ind,currSize_-1,
true);
249 L_.reshape(currSize_-1,currSize_-1);
253 unsigned zsize = currSize_-dependent_;
254 for(
unsigned m=ind; m<zsize; m++ ){
262 base_.erase(base_.begin()+ind);
271 Real ghNorm = GiGj(base_[currSize_-dependent_],base_[currSize_-dependent_]);
274 for (
unsigned ii=0; ii<currSize_-dependent_; ++ii){
275 lhnrm += L_(currSize_-dependent_,ii)*L_(currSize_-dependent_,ii);
277 deltaLh_ = std::abs(ghNorm - lhnrm);
279 Real sgn1 = sgn(deltaLh_);
280 Real sgn2 = sgn(dletaLj);
281 Real signum = sgn1 * sgn2;
282 deltaLh_ = std::abs( ghNorm - lhNorm + deltaLh_ * deltaLh_);
285 if( std::sqrt(deltaLh_) > tol*kappa_*std::max(static_cast<Real>(1),ghNorm) ){
286 unsigned newind = currSize_-dependent_;
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];
297 deltaLh_ = std::sqrt(deltaLh_);
298 addSubgradToBase(newind,deltaLh_);
302 Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
303 ljNorm -= deltaLj_ * deltaLj_;
305 deltaLj_ = std::abs(gjNorm - ljNorm);
307 deltaLj_ = - std::sqrt( deltaLj_ );
309 deltaLj_ = std::sqrt( deltaLj_ );
312 Real gjTgh = GiGj(base_[currSize_-1],base_[currSize_-2]);
314 for (
unsigned ii=0;ii<currSize_;ii++){
315 ljTlh += L_(currSize_-1,ii)*L_(currSize_-2,ii);
317 deltaLj_ = (gjTgh - ljTlh) / deltaLh_;
319 L_(currSize_-1,currSize_-2) = deltaLj_;
325 Real gjNorm = GiGj(base_[currSize_-1],base_[currSize_-1]);
328 for (
unsigned ii=0; ii<currSize_; ++ii) {
329 ljnrm += L_(currSize_-1,ii)*L_(currSize_-1,ii);
331 deltaLj_ = std::abs(gjNorm - ljnrm);
333 deltaLj_ = std::abs( gjNorm - ljNorm + deltaLj_ * deltaLj_);
336 if( std::sqrt(deltaLj_) > tol*kappa_*std::max(static_cast<Real>(1),gjNorm) ){
337 unsigned newind = currSize_-1;
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];
348 deltaLj_ = std::sqrt(deltaLj_);
349 addSubgradToBase(newind,deltaLj_);
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);
355 deltaLh_ /= deltaLj_;
358 deltaLh_ = - std::sqrt( deltaLh_ );
361 deltaLh_ = std::sqrt ( deltaLh_ );
364 L_(currSize_-1,currSize_-2) = deltaLh_;
370 template<
typename Real>
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;
381 lapack_.TRTRS(
'L', tran,
'N', size, 1, L.values(), L.stride(), v.values(), v.stride(), &info );
385 template<
typename Real>
388 for(
int i=0;i<v.numRows();i++){
396 template<
typename Real>
398 const Real
zero(0), half(0.5), one(1);
399 Real z1z2(0), z1z1(0);
406 rho_ = ROL_INF<Real>();
411 L_(0,0) = std::sqrt(GiGj(0,0));
419 z1_.resize(1); z2_.resize(1);
420 z1_[0] = one/L_(0,0);
425 objval_ = ROL_INF<Real>();
426 minobjval_ = ROL_INF<Real>();
430 for (iter=0;iter<maxit;iter++){
433 switch( dependent_ ){
441 rho_ = ( one + z1z2/t )/z1z1;
442 tempv_ = z1_; tempv_.scale(rho_);
443 tempw1_ = z2_; tempw1_.scale(one/t);
445 solveSystem(currSize_,
'T',L_,tempv_);
455 LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-1,currSize_-1);
456 lh_.size(currSize_-1);
459 for(
unsigned i=0; i<currSize_-1; ++i){
460 Real tmp = L_(currSize_-1,i);
466 if (std::abs(lhz1_-one) <= tol*kappa_){
469 solveSystem(currSize_-1,
'T',LBprime,tempv_);
470 tempv_.resize(currSize_);
471 tempv_[currSize_-1] = -one;
479 Real tmp = ( one + z1z2 / t - rho_ * z1z1 )/( one - lhz1_ );
480 tempw1_ = z1_; tempw1_.scale(rho_);
481 tempw2_ = z2_; tempw2_.scale(one/t);
483 tempw2_ = lh_; tempw2_.scale(tmp);
485 solveSystem(currSize_-1,
'T',LBprime,tempw1_);
487 tempv_.resize(currSize_);
488 tempv_[currSize_-1] = tmp;
499 LA::Matrix<Real> LBprime( LA::Copy,L_,currSize_-2,currSize_-2 );
500 lj_.size(currSize_-2);
501 lh_.size(currSize_-2);
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);
512 if(std::abs(ljz1_-one) <= tol*kappa_){
515 solveSystem(currSize_-2,
'T',LBprime,tempv_);
516 tempv_.resize(currSize_);
517 tempv_[currSize_-2] =
zero;
518 tempv_[currSize_-1] = -one;
522 Real mu = ( one - lhz1_ )/( one - ljz1_ );
523 tempw1_ = lj_; tempw1_.scale(-mu);
525 solveSystem(currSize_-2,
'T',LBprime,tempw1_);
527 tempv_.resize(currSize_);
528 tempv_[currSize_-2] = -one;
529 tempv_[currSize_-1] = mu;
538 if (isFeasible(tempv_,tol*currSize_)){
541 for (
unsigned i=0; i<currSize_; ++i){
547 for (
unsigned i=0; i<currSize_; ++i){
555 if ( tempv_[currSize_-1] <
zero )
560 for(
unsigned kk=0; kk<currSize_; ++kk)
569 Real myeps = tol*currSize_;
570 Real step = ROL_INF<Real>();
571 for (
unsigned h=0; h<currSize_; ++h){
578 if (step <=
zero || step == ROL_INF<Real>()){
583 for (
unsigned i=0; i<currSize_; ++i)
590 bool deleted = optimal_;
591 for (
unsigned baseitem=0; baseitem<currSize_; ++baseitem){
597 if( base_[baseitem] == entering_ ){
598 taboo_.push_back(entering_);
604 deleteSubgradFromBase(baseitem,tol);
615 Real newobjval(0), Lin(0), Quad(0);
616 for (
unsigned i=0; i<currSize_; ++i){
619 if (rho_ == ROL_NINF<Real>()){
621 newobjval = -half*Quad;
625 newobjval = half*(rho_ + Lin/t);
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_);
641 if ( objval_ + std::max( tol*std::abs(objval_), ROL_EPSILON<Real>() ) <= minobjval_ ){
643 minobjval_ = objval_;
647 if ( (rho_ >= ROL_NINF<Real>()) && (objval_ <= ROL_NINF<Real>()) )
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>();
656 for (
unsigned bundleitem=0; bundleitem<Bundle_U<Real>::size(); ++bundleitem){
658 if ( std::find(taboo_.begin(),taboo_.end(),bundleitem) != taboo_.end() ){
662 if ( std::find(base_.begin(),base_.end(),bundleitem) == base_.end() ){
665 for (
unsigned j=0;j<currSize_;j++){
671 if (rho_ >= ROL_NINF<Real>()){
675 ro = ROL_NINF<Real>();
676 minobjval_ = ROL_INF<Real>();
677 objval_ = ROL_INF<Real>();
681 #if ( FIRST_VIOLATED )
682 entering_ = bundleitem;
686 if ( diff > olddiff ){
687 entering_ = bundleitem;
697 if (entering_ < maxind_){
700 base_.push_back(entering_);
702 unsigned zsize = currSize_ - dependent_;
706 for (
unsigned i=0; i<zsize; ++i){
707 lh_[i] = GiGj(entering_,base_[i]);
709 LA::Matrix<Real> LBprime( LA::Copy,L_,zsize,zsize);
710 solveSystem(zsize,
'N',LBprime,lh_);
711 for (
unsigned i=0; i<zsize; ++i){
712 lhz1_ += lh_[i]*z1_[i];
713 lhz2_ += lh_[i]*z2_[i];
716 Real nrm = lh_.dot(lh_);
717 Real delta = GiGj(entering_,entering_) - nrm;
727 L_.reshape(currSize_,currSize_);
728 zsize = currSize_ - dependent_;
729 for (
unsigned i=0; i<zsize-1; ++i){
730 L_(currSize_-1,i) = lh_[i];
733 Real deltaeps = tol*kappa_*std::max(one,lh_.dot(lh_));
734 if ( delta > deltaeps ){
736 unsigned ind = currSize_-1;
737 delta = std::sqrt(delta);
738 addSubgradToBase(ind,delta);
740 else if(delta < -deltaeps){
750 if( objval_ - std::max( tol*std::abs( objval_ ), ROL_EPSILON<Real>() ) > minobjval_ ){
764 template<
typename Real>
767 unsigned outermaxit = 20;
768 bool increase =
false, decrease =
false;
770 for (
unsigned it=0; it < outermaxit; ++it ){
771 iter += solveDual_TT(t,maxit,mytol);
772 if ( QPStatus_ == 1 ) {
775 else if ( QPStatus_ == -2 || QPStatus_ == -3 ) {
776 mytol /=
static_cast<Real
>(10);
780 mytol *=
static_cast<Real
>(10);
783 if ( (mytol > static_cast<Real>(1e-4))
784 || (mytol < static_cast<Real>(1e-16)) ){
787 if ( increase && decrease ) {
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
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)