MoochoPack_ReducedHessianSecantUpdateBFGSProjected_Strategy.cpp

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

Generated on Tue Jul 13 09:29:32 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7