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