MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_InitFinDiffReducedHessian_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 <math.h>
00030 
00031 #include <ostream>
00032 
00033 #include "MoochoPack_InitFinDiffReducedHessian_Step.hpp"
00034 #include "MoochoPack_moocho_algo_conversion.hpp"
00035 #include "IterationPack_print_algorithm_step.hpp"
00036 #include "NLPInterfacePack_NLPObjGrad.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00038 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00039 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00040 //#include "AbstractLinAlgPack_SpVectorClass.hpp"
00041 //#include "AbstractLinAlgPack/src/max_near_feas_step.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00043 #include "AbstractLinAlgPack_VectorMutable.hpp"
00044 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00045 #include "AbstractLinAlgPack_VectorOut.hpp"
00046 #include "Teuchos_dyn_cast.hpp"
00047 
00048 namespace {
00049 template< class T >
00050 inline
00051 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00052 } // end namespace
00053 
00054 namespace MoochoPack {
00055 
00056 InitFinDiffReducedHessian_Step::InitFinDiffReducedHessian_Step(
00057   EInitializationMethod   initialization_method
00058   ,value_type             max_cond
00059   ,value_type             min_diag
00060   ,value_type             step_scale
00061   )
00062   :initialization_method_(initialization_method)
00063   ,max_cond_(max_cond)
00064   ,min_diag_(min_diag)
00065   ,step_scale_(step_scale)
00066 {}
00067 
00068 bool InitFinDiffReducedHessian_Step::do_step(
00069   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00070   ,poss_type assoc_step_poss
00071   )
00072 {
00073   using Teuchos::dyn_cast;
00074   using LinAlgOpPack::Vt_S;
00075   using LinAlgOpPack::Vp_StV;
00076   using LinAlgOpPack::V_MtV;
00077   using AbstractLinAlgPack::max_near_feas_step;
00078 
00079   NLPAlgo          &algo  = rsqp_algo(_algo);
00080   NLPAlgoState         &s     = algo.rsqp_state();
00081   NLPObjGrad    &nlp   = dyn_cast<NLPObjGrad>(algo.nlp());
00082   
00083   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00084   std::ostream& out = algo.track().journal_out();
00085 
00086   // print step header.
00087   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00088     using IterationPack::print_algorithm_step;
00089     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00090   }
00091 
00092   // Get the iteration quantity container objects
00093   IterQuantityAccess<index_type>
00094     &num_basis_iq = s.num_basis();
00095 
00096   const bool new_basis     = num_basis_iq.updated_k(0);
00097   const int  k_last_offset = s.rHL().last_updated();
00098 
00099   // If the basis has changed or there is no previous matrix to use
00100   // then reinitialize.
00101 
00102   if( new_basis || k_last_offset == IterQuantity::NONE_UPDATED ) {
00103 
00104     // Compute a finite difference along the null space of the
00105     // constraints
00106 
00107     IterQuantityAccess<VectorMutable>
00108       &x_iq     = s.x(),
00109       &rGf_iq   = s.rGf();
00110     IterQuantityAccess<MatrixOp>
00111       &Z_iq     = s.Z();
00112     IterQuantityAccess<MatrixSymOp>
00113       &rHL_iq   = s.rHL();
00114 
00115     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00116       out << "\nReinitializing the reduced Hessain using a finite difference\n";
00117     }
00118 
00119     MatrixSymInitDiag &rHL_diag = dyn_cast<MatrixSymInitDiag>(rHL_iq.set_k(0));
00120     const MatrixOp    &Z_k      = Z_iq.get_k(0);
00121     const Vector    &x_k      = x_iq.get_k(0);
00122     const Vector    &rGf_k    = rGf_iq.get_k(0);
00123 
00124     // one vector
00125     VectorSpace::vec_mut_ptr_t  e = Z_k.space_rows().create_member(1.0);
00126 
00127     // Ze
00128     VectorSpace::vec_mut_ptr_t Ze = x_k.space().create_member();
00129     V_MtV( Ze.get(), Z_k, BLAS_Cpp::no_trans, *e );
00130     
00131     // This does not have to be an accurate finite difference so lets just
00132     // take step_scale/||Ze|| as the step size unless it is out of bounds
00133     // If we assume that our variables are scaled near
00134     // one (step_scale == 1?) then this will give us an appreciable step.  Beside we
00135     // should not be near the solution so the reduced gradient should not
00136     // be near zero.
00137 
00138     const value_type nrm_Ze = Ze->norm_inf();
00139     value_type u = step_scale() / nrm_Ze;
00140 
00141     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00142       out << "\n||Ze||inf = " << nrm_Ze << std::endl;
00143     }
00144 
00145     if( (int)olevel >= (int)PRINT_VECTORS ) {
00146       out << "\nZe =\n" << *Ze;
00147     }
00148 
00149     if( algo.nlp().num_bounded_x() ) {
00150 
00151       // Find the maximum step u
00152       // we can take along x_fd = x_k + u*Ze
00153       // that don't violate variable bounds by too much.
00154       // If a positive step can't be found then this may be a negative step.
00155       
00156       std::pair<value_type,value_type>
00157         u_steps = max_near_feas_step(
00158           x_k, *Ze
00159           ,nlp.xl(), nlp.xu()
00160           ,nlp.max_var_bounds_viol()
00161           );
00162       
00163       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00164         out << "\nMaximum steps ( possitive, negative ) in bounds u = ("
00165           << u_steps.first << "," << u_steps.second << ")\n";
00166       }
00167 
00168       if( u_steps.first < u )
00169         u = u_steps.first;
00170       if( ::fabs(u_steps.second) < u )
00171         u = u_steps.second;
00172     }
00173 
00174     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00175       out << "\nFinite difference step length u = " << u << std::endl;
00176     }
00177 
00178     // Take the finite difference from x_k to x_fd = x_k + u * Ze
00179     //
00180     // rGf_fd = ( Z_k'*Gf(x_k + u*Ze) - rGf_k ) / u
00181     //
00182 
00183     VectorSpace::vec_mut_ptr_t x_fd = x_k.space().create_member();
00184     Vp_StV( x_fd.get(), u, *Ze );
00185 
00186     // Gf_fd = Gf(x_fd)
00187     VectorSpace::vec_mut_ptr_t Gf_fd = x_k.space().create_member();
00188     nlp.unset_quantities();
00189     nlp.set_Gf( Gf_fd.get() );
00190     nlp.calc_Gf( *x_fd );
00191 
00192     if( (int)olevel >= (int)PRINT_VECTORS ) {
00193       out << "\nGf_fd =\n" << *Gf_fd;
00194     }
00195 
00196     // rGf_fd = Z'*Gf_fd
00197     VectorSpace::vec_mut_ptr_t rGf_fd = Z_k.space_rows().create_member();
00198     V_MtV( rGf_fd.get(), Z_k, BLAS_Cpp::trans, *Gf_fd );
00199 
00200     // rGf_fd = rGf_fd - rGf_k
00201     Vp_StV( rGf_fd.get(), -1.0, rGf_k );
00202 
00203     // rGf_fd = rGf_fd / u
00204     Vt_S( rGf_fd.get(), 1.0 / u );
00205 
00206     const value_type
00207       nrm_rGf_fd = rGf_fd->norm_inf();
00208 
00209     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00210       out << "\n||(rGf_fd - rGf_k)/u||inf = " << nrm_rGf_fd << std::endl;
00211     }
00212     if( (int)olevel >= (int)PRINT_VECTORS ) {
00213       out << "\n(rGf_fd - rGf_k)/u =\n" << *rGf_fd;
00214     }
00215 
00216     if( nrm_rGf_fd <= min_diag() ) {
00217       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00218         out << "\n||(rGf_fd - rGf_k)/u||inf = " << nrm_rGf_fd
00219             << " < min_diag = " << min_diag() << std::endl
00220           << "\nScale by min_diag ... \n";
00221       }
00222       rHL_diag.init_identity(Z_k.space_rows(),min_diag());
00223     }
00224     else {
00225       switch( initialization_method() ) {
00226         case SCALE_IDENTITY: {
00227           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00228             out << "\nScale the identity matrix by ||(rGf_fd - rGf_k)/u||inf ... \n";
00229           }
00230           rHL_diag.init_identity(Z_k.space_rows(),nrm_rGf_fd);
00231           break;
00232         }
00233         case SCALE_DIAGONAL:
00234         case SCALE_DIAGONAL_ABS:
00235         {
00236           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00237             out << "\nScale diagonal by modified finite difference ... \n";
00238           }
00239           // In order to keep the initial reduced Hessian well conditioned we
00240           // will not let any diagonal element drop below
00241           // ||rGf_fd||inf / max_cond
00242         
00243           const value_type
00244             min_ele = my_max( nrm_rGf_fd / max_cond(), min_diag() );
00245 
00246           if( initialization_method() == SCALE_DIAGONAL )
00247             AbstractLinAlgPack::max_vec_scalar( min_ele, rGf_fd.get() );
00248           else
00249             AbstractLinAlgPack::max_abs_vec_scalar( min_ele, rGf_fd.get() );
00250 
00251           if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00252             out << "\n||diag||inf = " << rGf_fd->norm_inf() << std::endl;
00253           }
00254           if( (int)olevel >= (int)PRINT_VECTORS ) {
00255             out << "\ndiag =\n" << *rGf_fd;
00256           }
00257           rHL_diag.init_diagonal(*rGf_fd);
00258           break;
00259         }
00260         default:
00261           TEST_FOR_EXCEPT(true);  // only local programming error?
00262       }
00263     }
00264     nlp.unset_quantities();
00265 
00266     quasi_newton_stats_(s).set_k(0).set_updated_stats(
00267       QuasiNewtonStats::REINITIALIZED );
00268 
00269     if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00270       out << "\nrHL_k =\n" << rHL_iq.get_k(0);
00271     }
00272 
00273   }
00274 
00275   return true;
00276 }
00277 
00278 void InitFinDiffReducedHessian_Step::print_step(
00279   const Algorithm& algo
00280   ,poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00281   ,std::ostream& out, const std::string& L
00282   ) const
00283 {
00284   out
00285     << L << "*** Initialize the reduced Hessian using a single finite difference.\n"
00286     << L << "*** Where the nlp must support the NLPObjGrad interface and\n"
00287     << L << "*** rHL_k must support the MatrixSymInitDiag interface or exceptions\n"
00288     << L << "*** will be thrown.\n"
00289     << L << "default: num_basis_remembered = NO_BASIS_UPDATED_YET\n"
00290     << L << "         initialization_method = SCALE_DIAGONAL\n"
00291     << L << "         max_cond              = " << max_cond() << std::endl
00292     << L << "         min_diag              = " << min_diag() << std::endl
00293     << L << "         step_scale            = " << step_scale() << std::endl
00294     << L << "if num_basis_k is updated then\n"
00295     << L << "    new_basis = true\n"
00296     << L << "else\n"
00297     << L << "    new_basis = false\n"
00298     << L << "end\n"
00299     << L << "if new_basis == true or no past rHL as been updated then\n"
00300     << L << "    *** Reinitialize the reduced Hessian using finite differences\n"
00301     << L << "    Ze = Z * e\n"
00302     << L << "    u = step_scale / norm(Ze,inf)\n"
00303     << L << "    if there are bounds on the problem then\n"
00304     << L << "        Find the largest (in magnitude) positive (u_pos) and\n"
00305     << L << "        negative (u_neg) step u where the slightly relaxed variable\n"
00306     << L << "        bounds:\n"
00307     << L << "            xl - delta <= x_k + u * Ze <= xu + delta\n"
00308     << L << "        are strictly satisfied (where delta = max_var_bounds_viol).\n"
00309     << L << "        if u_pos < u then\n"
00310     << L << "            u = u_pos\n"
00311     << L << "        end\n"
00312     << L << "        if abs(u_neg) < u then\n"
00313     << L << "            u = u_neg\n"
00314     << L << "        end\n"
00315     << L << "    end\n"
00316     << L << "    x_fd = x_k + u * Ze\n"
00317     << L << "    rGf_fd = ( Z_k' * Gf(x_fd) - rGf_k ) / u\n"
00318     << L << "    if norm(rGf_fd,inf) <= min_diag then\n"
00319     << L << "        rHL_k = min_diag * eye(n-r)\n"
00320     << L << "    else\n"
00321     << L << "        if initialization_method == SCALE_IDENTITY then\n"
00322     << L << "            rHL_k = norm(rGf_fd,inf) * eye(n-r)\n"
00323     << L << "        else if initialization_method == SCALE_DIAGONAL or SCALE_DIAGONAL_ABS then\n"
00324     << L << "            *** Make sure that all of the diagonal elements are\n"
00325     << L << "            *** positive and that the smallest element is\n"
00326     << L << "            *** no smaller than norm(rGf_fd,inf) / max_cond\n"
00327     << L << "            *** So that rHL_k will be positive definite an\n"
00328     << L << "            *** well conditioned\n"
00329     << L << "            min_ele = max( norm(rGf_fd,inf)/max_cond, min_diag )\n"
00330     << L << "            if initialization_method == SCALE_DIAGONAL then\n"
00331     << L << "                for i = 1 ... n-r\n"
00332     << L << "                   diag(i) = max( rGf_fd(i), min_ele )\n"
00333     << L << "                end\n"
00334     << L << "            else *** SCALE_DIAGONAL_ABS\n"
00335     << L << "                for i = 1 ... n-r\n"
00336     << L << "                   diag(i) = max( abs(rGf_fd(i)), min_ele )\n"
00337     << L << "                end\n"
00338     << L << "            end\n"
00339     << L << "            rHL_k = diag(diag)\n"
00340     << L << "        end\n"
00341     << L << "    end\n"
00342     << L << "end\n";
00343 }
00344 
00345 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends