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