MoochoPack_ReducedHessianSecantUpdateBFGSProjected_Strategy.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "MoochoPack_ReducedHessianSecantUpdateBFGSProjected_Strategy.hpp"
00030 #include "MoochoPack_PBFGS_helpers.hpp"
00031 #include "MoochoPack_NLPAlgo.hpp"
00032 #include "MoochoPack_NLPAlgoState.hpp"
00033 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixSymAddDelUpdateable.hpp"
00034 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
00035 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00036 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00037 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOpOut.hpp"
00038 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00039 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00040 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSliceOut.hpp"
00041 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00042 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Midynamic_cast_verbose.h"
00044 #include "MiWorkspacePack.h"
00045 
00046 namespace LinAlgOpPack {
00047   using AbstractLinAlgPack::Vp_StMtV;
00048 }
00049 
00050 namespace MoochoPack {
00051 
00052 ReducedHessianSecantUpdateBFGSProjected_Strategy::ReducedHessianSecantUpdateBFGSProjected_Strategy(
00053   const bfgs_update_ptr_t&      bfgs_update
00054   ,value_type                   act_set_frac_proj_start
00055   ,value_type                   project_error_tol
00056   ,value_type                   super_basic_mult_drop_tol
00057   )
00058   : bfgs_update_(bfgs_update)
00059   , act_set_frac_proj_start_(act_set_frac_proj_start)
00060   , project_error_tol_(project_error_tol)
00061   , super_basic_mult_drop_tol_(super_basic_mult_drop_tol)
00062 {}
00063 
00064 bool ReducedHessianSecantUpdateBFGSProjected_Strategy::perform_update(
00065   DVectorSlice* s_bfgs, DVectorSlice* y_bfgs, bool first_update
00066   ,std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
00067   ,MatrixOp *rHL_k
00068   )
00069 {
00070   using std::setw;
00071   using std::endl;
00072   using std::right;
00073   using Teuchos::dyn_cast;
00074   using DenseLinAlgPack::dot;
00075   using LinAlgOpPack::V_MtV;
00076   using AbstractLinAlgPack::norm_inf;
00077   typedef ConstrainedOptPack::MatrixHessianSuperBasic MHSB_t;
00078   using Teuchos::Workspace;
00079   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00080 
00081   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00082     out << "\n*** (PBFGS) Projected BFGS ...\n";
00083   }
00084 
00085 #ifdef _WINDOWS
00086   MHSB_t &rHL_super = dynamic_cast<MHSB_t&>(*rHL_k);
00087 #else
00088   MHSB_t &rHL_super = dyn_cast<MHSB_t>(*rHL_k);
00089 #endif
00090 
00091   const size_type
00092     n    = algo->nlp().n(),
00093     r    = algo->nlp().r(),
00094     n_pz = n-r;
00095   const GenPermMatrixSlice
00096     &Q_R = rHL_super.Q_R(),
00097     &Q_X = rHL_super.Q_X();
00098 
00099   bool do_projected_rHL_RR = false;
00100 
00101   // Check that the current update is sufficiently p.d. before we do anything
00102   const value_type
00103     sTy = dot(*s_bfgs,*y_bfgs),
00104     yTy = dot(*y_bfgs,*y_bfgs);
00105   if( !ConstrainedOptPack::BFGS_sTy_suff_p_d(
00106     *s_bfgs,*y_bfgs,&sTy
00107     ,  int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL )
00108     && !bfgs_update().use_dampening()
00109     )
00110   {
00111     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00112       out << "\nWarning!  use_damening == false so there is no way we can perform any kind BFGS update (projected or not) so we will skip it!\n";
00113     }
00114     quasi_newton_stats_(*s).set_k(0).set_updated_stats(
00115       QuasiNewtonStats::INDEF_SKIPED );
00116     return true;
00117   }
00118 
00119   // Get matrix scaling
00120   value_type
00121     rHL_XX_scale = sTy > 0.0 ? yTy/sTy : 1.0;
00122 
00123   // 
00124   // Initialize or adjust the active set before the BFGS update
00125   //
00126   if( !s->nu().updated_k(-1) ) {
00127     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00128       out << "\nWarning!  nu_k(-1) has not been updated.  No adjustment to the set of superbasic variables is possible!\n";
00129     }
00130   }
00131   else if( Q_R.is_identity() ) {
00132     // Determine when to start adding and removing rows/cols form rHL_RR
00133     if( act_set_stats_(*s).updated_k(-1) ) {
00134       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00135         out << "\nDetermining if projected BFGS updating of superbasics should begin ...\n";
00136       }
00137       // Determine if the active set has calmed down enough
00138       const SpVector
00139         &nu_km1 = s->nu().get_k(-1);
00140       const SpVectorSlice
00141         nu_indep = nu_km1(s->var_indep());
00142       do_projected_rHL_RR = PBFGSPack::act_set_calmed_down(
00143         act_set_stats_(*s).get_k(-1)
00144         ,act_set_frac_proj_start()
00145         ,olevel,out
00146         );
00147       if( do_projected_rHL_RR ) {
00148         //
00149         // Determine the set of initially fixed and free independent variables.
00150         //
00151         typedef Workspace<size_type>                              i_x_t;
00152         typedef Workspace<ConstrainedOptPack::EBounds>   bnd_t;
00153         i_x_t   i_x_free(wss,n_pz);
00154         i_x_t   i_x_fixed(wss,n_pz);
00155         bnd_t   bnd_fixed(wss,n_pz);
00156         i_x_t   l_x_fixed_sorted(wss,n_pz);
00157         size_type n_pz_X = 0, n_pz_R = 0;
00158         value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
00159         // Get the elements in i_x_free[] for variables that are definitely free
00160         // and initialize s_R'*y_R
00161         PBFGSPack::init_i_x_free_sRTsR_sRTyR(
00162           nu_indep, *s_bfgs, *y_bfgs
00163           , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR );  // We don't really want sRTRBBsR here
00164         // rHL_XX = some_scaling * I
00165         if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00166           out << "\nScaling for diagonal rHL_XX = rHL_XX_scale*I, rHL_XX_scale = " << rHL_XX_scale << std::endl;
00167         }
00168         Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
00169         DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
00170         rHL_XX_diag = rHL_XX_scale;
00171         // s_R'*B_RR_*s_R
00172         Workspace<value_type> Q_R_Q_RT_s_ws(wss,n_pz);
00173         DVectorSlice Q_R_Q_RT_s(&Q_R_Q_RT_s_ws[0],Q_R_Q_RT_s_ws.size());
00174         Q_R_Q_RT_s = 0.0;
00175         {for( size_type k = 0; k < n_pz_R; ++k ) {
00176           const size_type i = i_x_free[k];
00177           Q_R_Q_RT_s(i) = (*s_bfgs)(i);
00178         }}
00179         sRTBRRsR = AbstractLinAlgPack::transVtMtV( Q_R_Q_RT_s, *rHL_k, BLAS_Cpp::no_trans,  Q_R_Q_RT_s );
00180         // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
00181         // and initialize s_X'*B_XX*s_X and s_X*y_X
00182         PBFGSPack::sort_fixed_max_cond_viol(
00183           nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
00184           ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
00185         // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
00186         PBFGSPack::choose_fixed_free(
00187           project_error_tol(),super_basic_mult_drop_tol(),nu_indep
00188           ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
00189           ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
00190           ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
00191         //
00192         // Delete rows/cols from rHL_RR for fixed variables
00193         //
00194 #ifdef _WINDOWS
00195         MatrixSymAddDelUpdateable
00196           &rHL_RR = dynamic_cast<MatrixSymAddDelUpdateable&>(
00197             const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
00198             );
00199 #else
00200         MatrixSymAddDelUpdateable
00201           &rHL_RR = dyn_cast<MatrixSymAddDelUpdateable>(
00202             const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
00203             );
00204 #endif
00205         if( n_pz_R < n_pz ) {
00206           if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00207             out << "\nDeleting n_pz_X = " << n_pz_X << " rows/columns from rHL_RR for fixed independent variables...\n";
00208           }
00209           {for( size_type k = n_pz_X; k > 0; --k ) { // Delete from the largest to the smallest index (cheaper!)
00210             rHL_RR.delete_update( i_x_fixed[k-1], false );
00211           }}
00212           assert( rHL_RR.rows() == n_pz_R );
00213           if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00214             out << "\nrHL_RR after rows/columns where removed =\n" << *rHL_super.B_RR_ptr();
00215           }
00216           // Initialize rHL_XX = rHL_XX_scale*I
00217 #ifdef _WINDOWS
00218           MatrixSymInitDiag
00219             &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
00220               const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
00221               );
00222 #else
00223           MatrixSymInitDiag
00224             &rHL_XX = dyn_cast<MatrixSymInitDiag>(
00225               const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
00226               );
00227 #endif
00228           rHL_XX.init_identity(n_pz_X,rHL_XX_scale);
00229           // Reinitialize rHL for new active set
00230           rHL_super.initialize(
00231             n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
00232             ,rHL_super.B_RR_ptr(),NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
00233             );
00234           if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00235             out << "\nFull rHL after reinitialization but before BFGS update:\n"
00236               << "\nrHL =\n" << *rHL_k
00237               << "\nQ_R =\n" << rHL_super.Q_R()
00238               << "\nQ_X =\n" << rHL_super.Q_X();
00239           }
00240         }
00241         else {
00242           do_projected_rHL_RR = false;
00243           if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00244             out << "\nWith n_pz_X = " << n_pz_X << ", there where no variables to drop from superbasis!\n";
00245           }
00246         }
00247       }
00248     }
00249   }
00250   else {
00251     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00252       out << "\nAdjust the set of superbasic variables and the projected reduced Hessian rHL_RR ...\n";
00253     }
00254     //
00255     // Modify rHL_RR by adding and dropping rows/cols for freeded and fixed variables
00256     //
00257     const SpVectorSlice
00258       nu_indep = s->nu().get_k(-1)(s->var_indep());
00259     //
00260     // Determine new Q_R and Q_X
00261     //
00262     typedef Workspace<size_type>                              i_x_t;
00263     typedef Workspace<ConstrainedOptPack::EBounds>   bnd_t;
00264     i_x_t   i_x_free(wss,n_pz);
00265     i_x_t   i_x_fixed(wss,n_pz);
00266     bnd_t   bnd_fixed(wss,n_pz);
00267     i_x_t   l_x_fixed_sorted(wss,n_pz);
00268     size_type n_pz_X = 0, n_pz_R = 0;
00269     value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
00270     // Get the elements in i_x_free[] for variables that are definitely free
00271     // and initialize s_R'*y_R.  This will be the starting point for the new Q_R.
00272     PBFGSPack::init_i_x_free_sRTsR_sRTyR(
00273       nu_indep, *s_bfgs, *y_bfgs
00274       , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR );  // We don't really want sRTBRRsR here
00275     // Initialize rHL_XX_diag = some_scaling * I as though all of the currently fixed variables
00276     // will be left out of Q_R.  Some of these variables might already be in Q_R and B_RR
00277     // and may still be in Q_R and B_RR after we are finished adjusting the sets Q_R and Q_X
00278     // and we don't want to delete these rows/cols in B_RR just yet!
00279     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00280       out << "\nScaling for diagonal rHL_XX = rHL_XX_scale*I, rHL_XX_scale = " << rHL_XX_scale << std::endl;
00281     }
00282     Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
00283     DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
00284     rHL_XX_diag = rHL_XX_scale;
00285     // Initialize rHL_XX = rHL_XX_scale * I so that those variables in the current Q_R
00286     // not in the estimate i_x_free[] will have their proper value when s_R'*B_RR*s_R computed
00287     // for the estimate i_x_free[].  This is needed to change the behavior of *rHL_k which
00288     // is used below to compute s_R'*B_RR*s_R
00289 #ifdef _WINDOWS
00290     MatrixSymInitDiag
00291       &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
00292         const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
00293         );
00294 #else
00295     MatrixSymInitDiag
00296       &rHL_XX = dyn_cast<MatrixSymInitDiag>(
00297         const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr())
00298         );
00299 #endif
00300     rHL_XX.init_identity(rHL_XX.rows(),rHL_XX_scale); // Don't change its size yet!
00301     // s_R'*B_RR_*s_R
00302     // This will only include those terms for the variable actually free.
00303     Workspace<value_type> Q_R_Q_RT_s_ws(wss,n_pz);
00304     DVectorSlice Q_R_Q_RT_s(&Q_R_Q_RT_s_ws[0],Q_R_Q_RT_s_ws.size());
00305     Q_R_Q_RT_s = 0.0;
00306     {for( size_type k = 0; k < n_pz_R; ++k ) {
00307       const size_type i = i_x_free[k];
00308       Q_R_Q_RT_s(i) = (*s_bfgs)(i);
00309     }}
00310     sRTBRRsR = AbstractLinAlgPack::transVtMtV( Q_R_Q_RT_s, *rHL_k, BLAS_Cpp::no_trans,  Q_R_Q_RT_s );
00311     // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
00312     PBFGSPack::sort_fixed_max_cond_viol(
00313       nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
00314       ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
00315     // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
00316     PBFGSPack::choose_fixed_free(
00317       project_error_tol(),super_basic_mult_drop_tol(),nu_indep
00318       ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
00319       ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
00320       ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
00321     // Get the changes to the set of superbasic variables
00322     size_type num_free_to_fixed = 0, num_fixed_to_free = 0;
00323     i_x_t  i_x_free_to_fixed(wss,Q_R.cols());
00324     i_x_t  i_x_fixed_to_free(wss,Q_X.cols());
00325     i_x_t  i_x_free_still(wss,Q_R.cols());             // Will be set to those indices still in Q_R
00326     std::fill_n( &i_x_free_still[0], Q_R.cols(), 0 );  // in the same order as in Q_R and B_RR
00327     {
00328       GenPermMatrixSlice::const_iterator
00329         Q_R_begin     = Q_R.begin(),
00330         Q_R_itr       = Q_R_begin,
00331         Q_R_end       = Q_R.end();
00332       const size_type
00333         *i_x_free_itr = &i_x_free[0],
00334         *i_x_free_end = i_x_free_itr + n_pz_R;
00335       for( size_type i = 1; i <= n_pz; ++i ) {
00336         if( Q_R_itr == Q_R_end && i_x_free_itr == i_x_free_end ) {
00337           break; // The rest of these variables were and still are not superbasic
00338         }
00339         else if( i_x_free_itr == i_x_free_end ) {
00340           // A variable that was in the superbasis now is not
00341           i_x_free_to_fixed[num_free_to_fixed] = Q_R_itr->row_i();
00342           num_free_to_fixed++;
00343           ++Q_R_itr;
00344         }
00345         else if( Q_R_itr == Q_R_end ) {
00346           // A variable that was not in the superbasis now is
00347           i_x_fixed_to_free[num_fixed_to_free] = *i_x_free_itr;
00348           num_fixed_to_free++;
00349           ++i_x_free_itr;
00350         }
00351         else {
00352           if( Q_R_itr->row_i() == *i_x_free_itr ) {
00353             // Both still superbasic
00354             i_x_free_still[Q_R_itr-Q_R_begin] = Q_R_itr->row_i();
00355             ++Q_R_itr;
00356             ++i_x_free_itr;
00357           }
00358           else if( Q_R_itr->row_i() < *i_x_free_itr ) {
00359             // A variable that was in the superbasis now is not
00360             i_x_free_to_fixed[num_free_to_fixed] = Q_R_itr->row_i();
00361             num_free_to_fixed++;
00362             ++Q_R_itr;
00363           }
00364           else {
00365             // A variable that was not in the superbasis now is
00366             i_x_fixed_to_free[num_fixed_to_free] = *i_x_free_itr;
00367             num_fixed_to_free++;
00368             ++i_x_free_itr;
00369           }
00370         }
00371       }
00372     }
00373     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00374       out << "\nThere will be " << num_fixed_to_free  << " independent variables added to the superbasis and rHL_RR";
00375       if( num_fixed_to_free && int(olevel) >= int(PRINT_ACTIVE_SET) ) {
00376         out << " and their indexes are:\n";
00377         for(size_type k = 0; k < num_fixed_to_free; ++k)
00378           out << " " << i_x_fixed_to_free[k];
00379         out << std::endl;
00380       }
00381       else {
00382         out << "\n";
00383       }
00384       out << "\nThere will be " << num_free_to_fixed  << " independent variables removed from the superbasis and rHL_RR";
00385       if( num_free_to_fixed && int(olevel) >= int(PRINT_ACTIVE_SET) ) {
00386         out << " and their indexes are:\n";
00387         for(size_type k = 0; k < num_free_to_fixed; ++k)
00388           out << " " << i_x_free_to_fixed[k];
00389         out << std::endl;
00390       }
00391       else {
00392         out << "\n";
00393       }
00394     }
00395     // Get reference to rHL_RR = B_RR
00396 #ifdef _WINDOWS
00397     MatrixSymAddDelUpdateable
00398       &rHL_RR = dynamic_cast<MatrixSymAddDelUpdateable&>(
00399         const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
00400         );
00401 #else
00402     MatrixSymAddDelUpdateable
00403       &rHL_RR = dyn_cast<MatrixSymAddDelUpdateable>(
00404         const_cast<MatrixSymWithOpFactorized&>(*rHL_super.B_RR_ptr())
00405         );
00406 #endif
00407     // Remove rows/cols from rHL_RR.
00408     if( num_free_to_fixed ) {
00409       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00410         out << "\nDeleting " << num_free_to_fixed << " rows/columns from rHL_RR ...\n";
00411       }
00412       {for( size_type k = i_x_free_still.size(); k > 0; --k ) { // Delete from the largest to the smallest index (cheaper!)
00413         if( !i_x_free_still[k-1] )
00414           rHL_RR.delete_update( k, false );
00415       }}
00416       assert( rHL_RR.rows() == n_pz_R - num_fixed_to_free );
00417       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00418         out << "\nrHL_RR after rows/columns where removed =\n" << *rHL_super.B_RR_ptr();
00419       }
00420     }
00421     // Add new rows/cols to rHL_RR.
00422     if( num_fixed_to_free ) {
00423       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00424         out << "\nAppending " << num_fixed_to_free << " rows/columns to rHL_RR ...\n";
00425       }
00426       {for( size_type k = 0; k < num_fixed_to_free; ++k ) {
00427         rHL_RR.augment_update( NULL, rHL_XX_scale, false );
00428       }}
00429       assert( rHL_RR.rows() == n_pz_R );
00430       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00431         out << "\nrHL_RR after rows/columns where appended =\n" << *rHL_super.B_RR_ptr();
00432       }
00433     }
00434     // Resort i_x_free[] to reflect the actual order of the indices in rHL_RR
00435     {
00436       size_type tmp_n_pz_R = 0;
00437       {for(size_type k = 0; k < i_x_free_still.size(); ++k) {
00438         if( i_x_free_still[k] ) {
00439           i_x_free[tmp_n_pz_R] = i_x_free_still[k];
00440           ++tmp_n_pz_R;
00441         }
00442       }}
00443       {for(size_type k = 0; k < num_fixed_to_free; ++k) {
00444         i_x_free[tmp_n_pz_R] = i_x_fixed_to_free[k];
00445         ++tmp_n_pz_R;
00446       }}
00447       assert( tmp_n_pz_R == n_pz_R );
00448     }
00449     // Initialize rHL_XX = rHL_XX_scale * I resized to the proper dimensions
00450     rHL_XX.init_identity(n_pz_X,rHL_XX_scale);
00451     // Reinitalize rHL for new active set
00452     rHL_super.initialize(
00453       n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
00454       ,rHL_super.B_RR_ptr(),NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
00455       );
00456     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00457       out << "\nFull rHL after reinitialization but before BFGS update:\n"
00458         << "\nrHL =\n" << *rHL_k
00459         << "\nQ_R =\n" << rHL_super.Q_R()
00460         << "\nQ_X =\n" << rHL_super.Q_X();
00461     }
00462     // Now we will do the PBFGS updating from now on!
00463     do_projected_rHL_RR = true;
00464   }
00465   //
00466   // Perform the BFGS update
00467   //
00468   if( do_projected_rHL_RR ) {
00469     // Perform BFGS update on smaller rHL_RR.
00470     // By the time we get here rHL_RR should be resize and ready to update
00471     const GenPermMatrixSlice
00472       &Q_R = rHL_super.Q_R(),
00473       &Q_X = rHL_super.Q_X();
00474     const size_type
00475       n_pz_R = Q_R.cols(),
00476       n_pz_X = Q_X.cols();
00477     assert( n_pz_R + n_pz_X == n_pz );
00478     // Get projected BFGS update vectors y_bfgs_R, s_bfgs_R
00479     Workspace<value_type>
00480       y_bfgs_R_ws(wss,Q_R.cols()),
00481       s_bfgs_R_ws(wss,Q_R.cols());
00482     DVectorSlice y_bfgs_R(&y_bfgs_R_ws[0],y_bfgs_R_ws.size());
00483     DVectorSlice s_bfgs_R(&s_bfgs_R_ws[0],s_bfgs_R_ws.size());
00484     V_MtV( &y_bfgs_R, Q_R, BLAS_Cpp::trans, *y_bfgs );  // y_bfgs_R = Q_R'*y_bfgs
00485     V_MtV( &s_bfgs_R, Q_R, BLAS_Cpp::trans, *s_bfgs );  // s_bfgs_R = Q_R'*s_bfgs
00486     // Update rHL_RR
00487     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00488       out << "\nPerform BFGS update on " << n_pz_R << " x " << n_pz_R << " projected reduced Hessian for the superbasic variables where B = rHL_RR...\n";
00489     }
00490     bfgs_update().perform_update(
00491       &s_bfgs_R(),&y_bfgs_R(),first_update,out,olevel,algo->algo_cntr().check_results()
00492       ,const_cast<MatrixOp*>(static_cast<const MatrixOp*>(rHL_super.B_RR_ptr().get()))
00493       ,&quasi_newton_stats_(*s).set_k(0)
00494       );
00495   }
00496   else {
00497     // Update the full reduced Hessain matrix (rHL = rHL_RR)
00498     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00499       out << "\nPerform BFGS update on the full reduced Hessian where B = rHL...\n";
00500     }
00501     bfgs_update().perform_update(
00502       s_bfgs,y_bfgs,first_update,out,olevel,algo->algo_cntr().check_results()
00503       ,const_cast<MatrixOp*>(static_cast<const MatrixOp*>(rHL_super.B_RR_ptr().get()))
00504       ,&quasi_newton_stats_(*s).set_k(0)
00505       );
00506   }
00507 
00508   return true;
00509 }
00510 
00511 void ReducedHessianSecantUpdateBFGSProjected_Strategy::print_step( std::ostream& out, const std::string& L ) const
00512 {
00513   out
00514     << L << "*** Perform BFGS update on only free independent (superbasic) variables.\n"
00515     << L << "if s'*y < sqrt(macheps) * ||s||2 * ||y||2 and use_dampening == false then"
00516     << L << "    Skip the update and exit this step!\n"
00517     << L << "end\n"
00518     << L << "rHL_XX_scale = max(y'*y/s'*y,1.0)\n"
00519     << L << "do_projected_rHL_RR = false\n"
00520     << L << "if nu_km1 is updated then\n"
00521     << L << "    if rHL_k.Q_R is the identity matrix then\n"
00522     << L << "        *** Determine if the active set has calmed down enough\n"
00523     << L << "        Given (num_active_indep,num_adds_indep,num_drops_indep) from act_set_stats_km1\n"
00524     << L << "        fact_same =\n"
00525     << L << "            ( num_adds_indep== NOT_KNOWN || num_drops_indep==NOT_KNOWN\n"
00526     << L << "            || num_active_indep==0\n"
00527     << L << "                ? 0.0\n"
00528     << L << "                : std::_MAX((num_active_indep-num_adds_indep-num_drops_indep)\n"
00529     << L << "                    /num_active_indep,0.0)\n"
00530     << L << "            )\n"
00531     << L << "        do_projected_rHL_RR\n"
00532     << L << "            = fact_same >= act_set_frac_proj_start && num_active_indep > 0\n"
00533     << L << "        if do_projected_rHL_RR == true then\n"
00534     << L << "            Determine the sets of superbasic variables given the mapping matrix\n"
00535     << L << "            Q = [ Q_R, Q_X ] where pz_R = Q_R'*pz <: R^n_pz_R are the superbasic variables and\n"
00536     << L << "            pz_X = Q_X'*pz <: R^n_pz_X are the nonbasic variables that only contain fixed\n"
00537     << L << "            variables in nu_km1(indep) where the following condidtions are satisfied:\n"
00538     << L << "                (s'*Q_X*rHL_XX_scale*Q_X'*s)/(s'*Q_R*Q_R'*rHL_k*Q_R*Q_R'*s) <= project_error_tol\n"
00539     << L << "                |s'*Q_X*Q_X'*y|/|s'*Q_R*Q_R'*s| <= project_error_tol\n"
00540     << L << "                |Q_X'*nu_km1(indep)|/||nu_km1(indep)||inf >= super_basic_mult_drop_tol\n"
00541     << L << "            if n_pz_R < n-r then\n"
00542     << L << "                Delete rows/columns of rHL_k to form rHL_RR = Q_R'*rHL_k*Q_R\n"
00543     << L << "                Define new rHL_k = [ Q_R, Q_X ] * [ rHL_RR, 0; 0; rHL_XX_scale*I ] [ Q_R'; Q_X ]\n"
00544     << L << "            else\n"
00545     << L << "                do_projected_rHL_RR = false\n"
00546     << L << "            end\n"
00547     << L << "        end\n"
00548     << L << "    else\n"
00549     << L << "        Determine the new Q_n = [ Q_R_n, Q_X_n ] that satisfies:\n"
00550     << L << "            (s'*Q_X_n*rHL_XX_scale*Q_X_n'*s)/(s'*Q_R_n*Q_R_n'*rHL_k*Q_R_n*Q_R_n'*s) <= project_error_tol\n"
00551     << L << "            |s'*Q_X_n*Q_X_n'*y|/|s'*Q_R_n*Q_R_n'*s| <= project_error_tol\n"
00552     << L << "            |Q_X_n'*nu_km1(indep)|/||nu_km1(indep)||inf >= super_basic_mult_drop_tol\n"
00553     << L << "        Remove rows/cols from rHL_k.rHL_RR for variables in rHL_k.Q_R that are not in Q_R_n.\n"
00554     << L << "        Add digonal entries equal to rHL_XX_scale to rHL_k.rHL_RR for variables in Q_R_n\n"
00555     << L << "        that are not in rHL_k.Q_R\n"
00556     << L << "        Define new rHL_k = [ Q_R_n, Q_X_n ] * [ rHL_k.rHL_RR, 0; 0; rHL_XX_scale*I ] [ Q_R_n'; Q_X_n ]\n"
00557     << L << "        do_projected_rHL_RR = true\n"
00558     << L << "    end\n"
00559     << L << "end\n"
00560     << L << "if do_projected_rHL_RR == true then\n"
00561     << L << "    Perform projected BFGS update (see below): (rHL_k.rHL_RR, Q_R_n'*s, Q_R_n'*y) -> rHL_k.rHL_RR\n"
00562     << L << "else\n"
00563     << L << "    Perform full BFGS update: (rHL_k, s, y) -> rHL_k\n"
00564     << L << "    begin BFGS update where B = rHL_k\n";
00565   bfgs_update().print_step( out, L + "        " );
00566   out
00567     << L << "    end BFGS update\n"
00568     << L << "else\n"
00569     ;
00570 }
00571 
00572 }  // end namespace MoochoPack

Generated on Thu Sep 18 12:35:18 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1