MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <ostream>
00043 
00044 #include "MoochoPack_ReducedHessianSecantUpdateStd_Step.hpp"
00045 #include "MoochoPack_moocho_algo_conversion.hpp"
00046 #include "MoochoPack_Exceptions.hpp"
00047 #include "IterationPack_print_algorithm_step.hpp"
00048 #include "AbstractLinAlgPack_VectorSpace.hpp"
00049 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00050 #include "AbstractLinAlgPack_VectorMutable.hpp"
00051 #include "AbstractLinAlgPack_VectorOut.hpp"
00052 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00053 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00054 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00056 #include "Teuchos_dyn_cast.hpp"
00057 
00058 MoochoPack::ReducedHessianSecantUpdateStd_Step::ReducedHessianSecantUpdateStd_Step(
00059   const secant_update_ptr_t&   secant_update
00060   )
00061   :secant_update_(secant_update)
00062   ,num_basis_(NO_BASIS_UPDATED_YET)
00063   ,iter_k_rHL_init_ident_(-1) // not a valid iteration?
00064 {}
00065 
00066 bool MoochoPack::ReducedHessianSecantUpdateStd_Step::do_step(
00067   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00068   ,poss_type assoc_step_poss
00069   )
00070 {
00071   using Teuchos::dyn_cast;
00072   using AbstractLinAlgPack::dot;
00073   using AbstractLinAlgPack::Vt_S;
00074   using AbstractLinAlgPack::Vp_StV;
00075   using LinAlgOpPack::Vp_V;
00076   using LinAlgOpPack::V_VpV;
00077   using LinAlgOpPack::V_VmV;
00078   using LinAlgOpPack::V_StV;
00079   using LinAlgOpPack::V_MtV;
00080 
00081   NLPAlgo       &algo = rsqp_algo(_algo);
00082   NLPAlgoState  &s    = algo.rsqp_state();
00083 
00084   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00085   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00086   std::ostream& out = algo.track().journal_out();
00087 
00088   // print step header.
00089   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00090     using IterationPack::print_algorithm_step;
00091     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00092   }
00093 
00094   bool return_val = true;
00095   
00096   // Get iteration quantities
00097   IterQuantityAccess<index_type>
00098     &num_basis_iq = s.num_basis();
00099   IterQuantityAccess<VectorMutable>
00100     &pz_iq  = s.pz(),
00101     &rGf_iq = s.rGf(),
00102     &w_iq   = s.w();
00103   IterQuantityAccess<MatrixOp>
00104     &Z_iq = s.Z();
00105   IterQuantityAccess<MatrixSymOp>
00106     &rHL_iq = s.rHL();
00107 
00108   // problem size
00109   const NLP &nlp = algo.nlp();
00110   const size_type
00111     //n    = nlp.n(),
00112     m    = nlp.m(),
00113     nind = m ? Z_iq.get_k(Z_iq.last_updated()).cols() : 0;
00114     //r    = m - nind;
00115 
00116   // See if a new basis has been selected
00117   bool new_basis = false;
00118   {
00119     const int last_updated_k = num_basis_iq.last_updated();
00120     if( last_updated_k != IterQuantity::NONE_UPDATED ) {
00121       const index_type num_basis_last = num_basis_iq.get_k(last_updated_k);
00122       if( num_basis_ == NO_BASIS_UPDATED_YET )
00123         num_basis_ = num_basis_last;
00124       else if( num_basis_ != num_basis_last )
00125         new_basis = true;
00126       num_basis_ = num_basis_last;
00127     }
00128   }
00129 
00130   // If rHL has already been updated for this iteration then just leave it.
00131   if( !rHL_iq.updated_k(0) ) {
00132 
00133     // If a new basis has been selected, reinitialize
00134     if( new_basis ) {
00135 
00136       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00137         out << "\nBasis changed.  Reinitializing rHL_k = eye(n-r) ...\n";
00138       }
00139       dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
00140         Z_iq.get_k(0).space_rows()
00141         );
00142       if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
00143         if( algo.algo_cntr().calc_matrix_norms() )
00144           out << "\n||rHL_k||inf = " << rHL_iq.get_k(0).calc_norm(MatrixOp::MAT_NORM_INF).value << std::endl;
00145       if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES )
00146         out << "\nrHL_k = \n" << rHL_iq.get_k(0);
00147       quasi_newton_stats_(s).set_k(0).set_updated_stats(
00148         QuasiNewtonStats::REINITIALIZED );
00149       iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
00150 
00151     }
00152     else {
00153       
00154       // Determine if rHL has been initialized and if we
00155       // can perform the update.  To perform the BFGS update
00156       // rHL_km1 and rGf_km1 must have been computed.
00157       if( rHL_iq.updated_k(-1) && rGf_iq.updated_k(-1) ) {
00158 
00159         // /////////////////////////////////////////////////////
00160         // Perform the Secant update
00161 
00162         if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
00163         {
00164           out
00165             << "\nPerforming Secant update ...\n";
00166         }
00167 
00168         const Vector
00169           &rGf_k   = rGf_iq.get_k(0),
00170           &rGf_km1 = rGf_iq.get_k(-1),
00171           &pz_km1  = pz_iq.get_k(-1);
00172         const value_type
00173           alpha_km1 = s.alpha().get_k(-1);
00174         VectorSpace::vec_mut_ptr_t
00175           y_bfgs = rGf_k.space().create_member(),
00176           s_bfgs = pz_km1.space().create_member();
00177 
00178         // /////////////////////////////////////////////////////
00179         // y_bfgs = rGf_k - rGf_km1 - alpha_km1 * w_km1
00180       
00181         // y_bfgs = rGf_k - rGf_km1 
00182         V_VmV( y_bfgs.get(), rGf_k, rGf_km1 );  
00183 
00184         if( w_iq.updated_k(-1) )
00185           // y_bfgs += - alpha_km1 * w_km1
00186           Vp_StV( y_bfgs.get(), - alpha_km1, w_iq.get_k(-1) );
00187 
00188         // /////////////////////////////////////////////////////
00189         // s_bfgs = alpha_km1 * pz_km1
00190         V_StV( s_bfgs.get(), alpha_km1, pz_iq.get_k(-1) );
00191 
00192         if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00193           out << "\n||y_bfgs||inf = " << y_bfgs->norm_inf() << std::endl;
00194           out << "\n||s_bfgs||inf = " << s_bfgs->norm_inf() << std::endl;
00195         }
00196 
00197         if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00198           out << "\ny_bfgs =\n" << *y_bfgs;
00199           out << "\ns_bfgs =\n" << *s_bfgs;
00200         }
00201 
00202         // Update from last
00203         MatrixSymOp
00204           &rHL_k   = rHL_iq.set_k(0,-1);
00205 
00206         // Perform the secant update
00207         if(!secant_update().perform_update(
00208              s_bfgs.get(), y_bfgs.get()
00209              ,iter_k_rHL_init_ident_ == s.k() - 1
00210              ,out, ns_olevel, &algo, &s, &rHL_k
00211              ))
00212         {
00213           return_val = false; // redirect control of algorithm!
00214         }
00215 
00216       }
00217       else {
00218         // We do not have the info to perform the update
00219 
00220         int k_last_offset = rHL_iq.last_updated();
00221         bool set_current = false;
00222         if( k_last_offset != IterQuantity::NONE_UPDATED && k_last_offset < 0 ) {
00223           const MatrixSymOp &rHL_k_last = rHL_iq.get_k(k_last_offset);
00224           const size_type nind_last = rHL_k_last.rows();
00225           if( nind_last != nind) {
00226             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00227               out
00228                 << "No new basis was selected.\n"
00229                 << "The previous matrix rHL_k(" << k_last_offset << ") was found but its dimmension\n"
00230                 << "rHL_k(" << k_last_offset << ").rows() = " << nind_last << " != n-r = " << nind << std::endl;
00231             }
00232           }
00233           else {
00234             if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00235               out
00236                 << "No new basis was selected so using previously updated...\n "
00237                 << "rHL_k = rHL_k(" << k_last_offset << ")\n";
00238             }
00239             rHL_iq.set_k(0) = rHL_k_last;
00240             quasi_newton_stats_(s).set_k(0).set_updated_stats(
00241               QuasiNewtonStats::SKIPED );
00242           
00243             if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00244               rHL_iq.get_k(0).output( out << "\nrHL_k = \n" );
00245             }
00246             set_current = true;
00247           }
00248         }
00249         if( !set_current ) {
00250           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00251             out
00252               << "\nInitializing rHL = eye(n-r) "
00253               << "(k = " << algo.state().k() << ")...\n";
00254           }
00255 
00256           // Now I will assume that since I can't perform the BFGS update and rHL has
00257           // not been set for this iteration yet, that it is up to me to initialize rHL_k = 0
00258           dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0)).init_identity(
00259             Z_iq.get_k(0).space_rows() );
00260           iter_k_rHL_init_ident_ = s.k(); // remember what iteration this was
00261           quasi_newton_stats_(s).set_k(0).set_updated_stats(
00262             QuasiNewtonStats::REINITIALIZED );
00263         
00264         }
00265       }
00266 
00267     }
00268     
00269     // Print rHL_k
00270     
00271     MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
00272 
00273     if( (int)ns_olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00274       if(algo.algo_cntr().calc_matrix_norms())
00275         out << "\n||rHL_k||inf    = " << rHL_iq.get_k(0).calc_norm(mat_nrm_inf).value << std::endl;
00276       if(algo.algo_cntr().calc_conditioning()) {
00277         const MatrixSymOpNonsing
00278           *rHL_ns_k = dynamic_cast<const MatrixSymOpNonsing*>(&rHL_iq.get_k(0));
00279         if(rHL_ns_k)
00280           out << "\ncond_inf(rHL_k) = " << rHL_ns_k->calc_cond_num(mat_nrm_inf).value << std::endl;
00281       }
00282     }
00283     if( (int)ns_olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00284       out << "\nrHL_k = \n" << rHL_iq.get_k(0);
00285     }
00286     
00287   }
00288   else {
00289     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00290       out << "\nThe matrix rHL_k has already been updated so leave it\n";
00291     }
00292   }
00293   
00294   return return_val;
00295 }
00296 
00297 void MoochoPack::ReducedHessianSecantUpdateStd_Step::print_step(
00298   const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00299   ,std::ostream& out, const std::string& L
00300   ) const
00301 {
00302   out
00303     << L << "*** Calculate the reduced hessian of the Lagrangian rHL = Z' * HL * Z\n"
00304     << L << "default:  num_basis_remembered = NO_BASIS_UPDATED_YET\n"
00305     << L << "          iter_k_rHL_init_ident = -1\n"
00306     << L << "if num_basis_remembered = NO_BASIS_UPDATED_YET then\n"
00307     << L << "    num_basis_remembered = num_basis\n"
00308     << L << "end\n"
00309     << L << "if num_basis_remembered != num_basis then\n"
00310     << L << "    num_basis_remembered = num_basis\n"
00311     << L << "    new_basis = true\n"
00312     << L << "end\n"
00313     << L << "if rHL_k is not updated then\n"
00314     << L << "    if new_basis == true then\n"
00315     << L << "        *** Transition rHL to the new basis by just starting over.\n"
00316     << L << "        rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
00317     << L << "        iter_k_rHL_init_ident = k\n"
00318     << L << "        goto next step\n"
00319     << L << "    end\n"
00320     << L << "    if rHL_km1 and rGf_km1 are updated then\n"
00321     << L << "        *** We should have the information to perform a BFGS update\n"
00322     << L << "        y = rGf_k - rGf_km1\n"
00323     << L << "        s = alpha_km1 * pz_km1\n"
00324     << L << "        if k - 1 == iter_k_rHL_init_ident then\n"
00325     << L << "            first_update = true\n"
00326     << L << "        else\n"
00327     << L << "            first_update = false\n"
00328     << L << "        end\n"
00329     << L << "        rHL_k = rHL_km1\n"
00330     << L << "        begin secant update\n"
00331     << L << "        (" << typeName(secant_update()) << ")\n"
00332     ;
00333   secant_update().print_step( out, L+"            " );
00334   out
00335     << L << "        end secant update\n"
00336     << L << "    else\n"
00337     << L << "       *** We have no information for which to preform a BFGS update.\n"
00338     << L << "       k_last_offset = last iteration rHL was updated for\n"
00339     << L << "       if k_last_offset does not exist then\n"
00340     << L << "            *** We are left with no choise but to initialize rHL\n"
00341     << L << "            rHL_k = eye(n-r) *** must support MatrixSymInitDiag interface\n"
00342     << L << "            iter_k_rHL_init_ident = k\n"
00343     << L << "        else\n"
00344     << L << "            *** No new basis has been selected so we may as well\n"
00345     << L << "            *** just use the last rHL that was updated\n"
00346     << L << "            rHL_k = rHL_k(k_last_offset)\n"
00347     << L << "        end\n"
00348     << L << "    end\n"
00349     << L << "end\n";
00350 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends