MoochoPack_ReducedHessianSecantUpdateStd_Step.cpp

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 <ostream>
00030 
00031 #include "MoochoPack_ReducedHessianSecantUpdateStd_Step.hpp"
00032 #include "MoochoPack_moocho_algo_conversion.hpp"
00033 #include "MoochoPack_Exceptions.hpp"
00034 #include "IterationPack_print_algorithm_step.hpp"
00035 #include "AbstractLinAlgPack_VectorSpace.hpp"
00036 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00037 #include "AbstractLinAlgPack_VectorMutable.hpp"
00038 #include "AbstractLinAlgPack_VectorOut.hpp"
00039 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00040 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00041 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "Teuchos_dyn_cast.hpp"
00044 
00045 MoochoPack::ReducedHessianSecantUpdateStd_Step::ReducedHessianSecantUpdateStd_Step(
00046   const secant_update_ptr_t&   secant_update
00047   )
00048   :secant_update_(secant_update)
00049   ,num_basis_(NO_BASIS_UPDATED_YET)
00050   ,iter_k_rHL_init_ident_(-1) // not a valid iteration?
00051 {}
00052 
00053 bool MoochoPack::ReducedHessianSecantUpdateStd_Step::do_step(
00054   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00055   ,poss_type assoc_step_poss
00056   )
00057 {
00058   using Teuchos::dyn_cast;
00059   using AbstractLinAlgPack::dot;
00060   using AbstractLinAlgPack::Vt_S;
00061   using AbstractLinAlgPack::Vp_StV;
00062   using LinAlgOpPack::Vp_V;
00063   using LinAlgOpPack::V_VpV;
00064   using LinAlgOpPack::V_VmV;
00065   using LinAlgOpPack::V_StV;
00066   using LinAlgOpPack::V_MtV;
00067 
00068   NLPAlgo       &algo = rsqp_algo(_algo);
00069   NLPAlgoState  &s    = algo.rsqp_state();
00070 
00071   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00072   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00073   std::ostream& out = algo.track().journal_out();
00074 
00075   // print step header.
00076   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00077     using IterationPack::print_algorithm_step;
00078     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00079   }
00080 
00081   bool return_val = true;
00082   
00083   // Get iteration quantities
00084   IterQuantityAccess<index_type>
00085     &num_basis_iq = s.num_basis();
00086   IterQuantityAccess<VectorMutable>
00087     &pz_iq  = s.pz(),
00088     &rGf_iq = s.rGf(),
00089     &w_iq   = s.w();
00090   IterQuantityAccess<MatrixOp>
00091     &Z_iq = s.Z();
00092   IterQuantityAccess<MatrixSymOp>
00093     &rHL_iq = s.rHL();
00094 
00095   // problem size
00096   const NLP &nlp = algo.nlp();
00097   const size_type
00098     //n    = nlp.n(),
00099     m    = nlp.m(),
00100     nind = m ? Z_iq.get_k(Z_iq.last_updated()).cols() : 0;
00101     //r    = m - nind;
00102 
00103   // See if a new basis has been selected
00104   bool new_basis = false;
00105   {
00106     const int last_updated_k = num_basis_iq.last_updated();
00107     if( last_updated_k != IterQuantity::NONE_UPDATED ) {
00108       const index_type num_basis_last = num_basis_iq.get_k(last_updated_k);
00109       if( num_basis_ == NO_BASIS_UPDATED_YET )
00110         num_basis_ = num_basis_last;
00111       else if( num_basis_ != num_basis_last )
00112         new_basis = true;
00113       num_basis_ = num_basis_last;
00114     }
00115   }
00116 
00117   // If rHL has already been updated for this iteration then just leave it.
00118   if( !rHL_iq.updated_k(0) ) {
00119 
00120     // If a new basis has been selected, reinitialize
00121     if( new_basis ) {
00122 
00123       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00124         out << "\nBasis changed.  Reinitializing rHL_k = eye(n-r) ...\n";
00125       }
00126       dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
00127         Z_iq.get_k(0).space_rows()
00128         );
00129       if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
00130         if( algo.algo_cntr().calc_matrix_norms() )
00131           out << "\n||rHL_k||inf = " << rHL_iq.get_k(0).calc_norm(MatrixOp::MAT_NORM_INF).value << std::endl;
00132       if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES )
00133         out << "\nrHL_k = \n" << rHL_iq.get_k(0);
00134       quasi_newton_stats_(s).set_k(0).set_updated_stats(
00135         QuasiNewtonStats::REINITIALIZED );
00136       iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
00137 
00138     }
00139     else {
00140       
00141       // Determine if rHL has been initialized and if we
00142       // can perform the update.  To perform the BFGS update
00143       // rHL_km1 and rGf_km1 must have been computed.
00144       if( rHL_iq.updated_k(-1) && rGf_iq.updated_k(-1) ) {
00145 
00146         // /////////////////////////////////////////////////////
00147         // Perform the Secant update
00148 
00149         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00150         {
00151           out
00152             << "\nPerforming Secant update ...\n";
00153         }
00154 
00155         const Vector
00156           &rGf_k   = rGf_iq.get_k(0),
00157           &rGf_km1 = rGf_iq.get_k(-1),
00158           &pz_km1  = pz_iq.get_k(-1);
00159         const value_type
00160           alpha_km1 = s.alpha().get_k(-1);
00161         VectorSpace::vec_mut_ptr_t
00162           y_bfgs = rGf_k.space().create_member(),
00163           s_bfgs = pz_km1.space().create_member();
00164 
00165         // /////////////////////////////////////////////////////
00166         // y_bfgs = rGf_k - rGf_km1 - alpha_km1 * w_km1
00167       
00168         // y_bfgs = rGf_k - rGf_km1 
00169         V_VmV( y_bfgs.get(), rGf_k, rGf_km1 );  
00170 
00171         if( w_iq.updated_k(-1) )
00172           // y_bfgs += - alpha_km1 * w_km1
00173           Vp_StV( y_bfgs.get(), - alpha_km1, w_iq.get_k(-1) );
00174 
00175         // /////////////////////////////////////////////////////
00176         // s_bfgs = alpha_km1 * pz_km1
00177         V_StV( s_bfgs.get(), alpha_km1, pz_iq.get_k(-1) );
00178 
00179         if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00180           out << "\n||y_bfgs||inf = " << y_bfgs->norm_inf() << std::endl;
00181           out << "\n||s_bfgs||inf = " << s_bfgs->norm_inf() << std::endl;
00182         }
00183 
00184         if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00185           out << "\ny_bfgs =\n" << *y_bfgs;
00186           out << "\ns_bfgs =\n" << *s_bfgs;
00187         }
00188 
00189         // Update from last
00190         MatrixSymOp
00191           &rHL_k   = rHL_iq.set_k(0,-1);
00192 
00193         // Perform the secant update
00194         if(!secant_update().perform_update(
00195              s_bfgs.get(), y_bfgs.get()
00196              ,iter_k_rHL_init_ident_ == s.k() - 1
00197              ,out, ns_olevel, &algo, &s, &rHL_k
00198              ))
00199         {
00200           return_val = false; // redirect control of algorithm!
00201         }
00202 
00203       }
00204       else {
00205         // We do not have the info to perform the update
00206 
00207         int k_last_offset = rHL_iq.last_updated();
00208         bool set_current = false;
00209         if( k_last_offset != IterQuantity::NONE_UPDATED && k_last_offset < 0 ) {
00210           const MatrixSymOp &rHL_k_last = rHL_iq.get_k(k_last_offset);
00211           const size_type nind_last = rHL_k_last.rows();
00212           if( nind_last != nind) {
00213             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00214               out
00215                 << "No new basis was selected.\n"
00216                 << "The previous matrix rHL_k(" << k_last_offset << ") was found but its dimmension\n"
00217                 << "rHL_k(" << k_last_offset << ").rows() = " << nind_last << " != n-r = " << nind << std::endl;
00218             }
00219           }
00220           else {
00221             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00222               out
00223                 << "No new basis was selected so using previously updated...\n "
00224                 << "rHL_k = rHL_k(" << k_last_offset << ")\n";
00225             }
00226             rHL_iq.set_k(0) = rHL_k_last;
00227             quasi_newton_stats_(s).set_k(0).set_updated_stats(
00228               QuasiNewtonStats::SKIPED );
00229           
00230             if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00231               rHL_iq.get_k(0).output( out << "\nrHL_k = \n" );
00232             }
00233             set_current = true;
00234           }
00235         }
00236         if( !set_current ) {
00237           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00238             out
00239               << "\nInitializing rHL = eye(n-r) "
00240               << "(k = " << algo.state().k() << ")...\n";
00241           }
00242 
00243           // Now I will assume that since I can't perform the BFGS update and rHL has
00244           // not been set for this iteration yet, that it is up to me to initialize rHL_k = 0
00245           dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
00246             Z_iq.get_k(0).space_rows() );
00247           iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
00248           quasi_newton_stats_(s).set_k(0).set_updated_stats(
00249             QuasiNewtonStats::REINITIALIZED );
00250         
00251         }
00252       }
00253 
00254     }
00255     
00256     // Print rHL_k
00257     
00258     MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
00259 
00260     if( (int)ns_olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00261       if(algo.algo_cntr().calc_matrix_norms())
00262         out << "\n||rHL_k||inf    = " << rHL_iq.get_k(0).calc_norm(mat_nrm_inf).value << std::endl;
00263       if(algo.algo_cntr().calc_conditioning()) {
00264         const MatrixSymOpNonsing
00265           *rHL_ns_k = dynamic_cast<const MatrixSymOpNonsing*>(&rHL_iq.get_k(0));
00266         if(rHL_ns_k)
00267           out << "\ncond_inf(rHL_k) = " << rHL_ns_k->calc_cond_num(mat_nrm_inf).value << std::endl;
00268       }
00269     }
00270     if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00271       out << "\nrHL_k = \n" << rHL_iq.get_k(0);
00272     }
00273     
00274   }
00275   else {
00276     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00277       out << "\nThe matrix rHL_k has already been updated so leave it\n";
00278     }
00279   }
00280   
00281   return return_val;
00282 }
00283 
00284 void MoochoPack::ReducedHessianSecantUpdateStd_Step::print_step(
00285   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00286   ,std::ostream& out, const std::string& L
00287   ) const
00288 {
00289   out
00290     << L << "*** Calculate the reduced hessian of the Lagrangian rHL = Z' * HL * Z\n"
00291     << L << "default:  num_basis_remembered = NO_BASIS_UPDATED_YET\n"
00292     << L << "          iter_k_rHL_init_ident = -1\n"
00293     << L << "if num_basis_remembered = NO_BASIS_UPDATED_YET then\n"
00294     << L << "    num_basis_remembered = num_basis\n"
00295     << L << "end\n"
00296     << L << "if num_basis_remembered != num_basis then\n"
00297     << L << "    num_basis_remembered = num_basis\n"
00298     << L << "    new_basis = true\n"
00299     << L << "end\n"
00300     << L << "if rHL_k is not updated then\n"
00301     << L << "    if new_basis == true then\n"
00302     << L << "        *** Transition rHL to the new basis by just starting over.\n"
00303     << L << "        rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
00304     << L << "        iter_k_rHL_init_ident = k\n"
00305     << L << "        goto next step\n"
00306     << L << "    end\n"
00307     << L << "    if rHL_km1 and rGf_km1 are updated then\n"
00308     << L << "        *** We should have the information to perform a BFGS update\n"
00309     << L << "        y = rGf_k - rGf_km1\n"
00310     << L << "        s = alpha_km1 * pz_km1\n"
00311     << L << "        if k - 1 == iter_k_rHL_init_ident then\n"
00312     << L << "            first_update = true\n"
00313     << L << "        else\n"
00314     << L << "            first_update = false\n"
00315     << L << "        end\n"
00316     << L << "        rHL_k = rHL_km1\n"
00317     << L << "        begin secant update\n"
00318     << L << "        (" << typeName(secant_update()) << ")\n"
00319     ;
00320   secant_update().print_step( out, L+"            " );
00321   out
00322     << L << "        end secant update\n"
00323     << L << "    else\n"
00324     << L << "       *** We have no information for which to preform a BFGS update.\n"
00325     << L << "       k_last_offset = last iteration rHL was updated for\n"
00326     << L << "       if k_last_offset does not exist then\n"
00327     << L << "            *** We are left with no choise but to initialize rHL\n"
00328     << L << "            rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
00329     << L << "            iter_k_rHL_init_ident = k\n"
00330     << L << "        else\n"
00331     << L << "            *** No new basis has been selected so we may as well\n"
00332     << L << "            *** just use the last rHL that was updated\n"
00333     << L << "            rHL_k = rHL_k(k_last_offset)\n"
00334     << L << "        end\n"
00335     << L << "    end\n"
00336     << L << "end\n";
00337 }

Generated on Wed May 12 21:51:20 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7