MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_MeritFunc_ModifiedL1LargerSteps_AddedStep.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 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include <ostream>
00032 #include <typeinfo>
00033 
00034 #include "MoochoPack_MeritFunc_ModifiedL1LargerSteps_AddedStep.hpp"
00035 #include "MoochoPack_moocho_algo_conversion.hpp"
00036 #include "IterationPack_print_algorithm_step.hpp"
00037 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
00038 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
00039 #include "ConstrainedOptPack/src/VectorWithNorms.h"
00040 #include "DenseLinAlgPack_DVectorOp.hpp"
00041 #include "DenseLinAlgPack_DVectorClass.hpp"
00042 #include "DenseLinAlgPack_DVectorOut.hpp"
00043 
00044 namespace {
00045 
00046 typedef MoochoPack::value_type value_type;
00047 inline value_type max(value_type v1, value_type v2)
00048 { return (v1 > v2) ? v1 : v2; }
00049 
00050 }
00051 
00052 namespace MoochoPack {
00053 
00054 MeritFunc_ModifiedL1LargerSteps_AddedStep::MeritFunc_ModifiedL1LargerSteps_AddedStep(
00055       const merit_func_ptr_t& merit_func
00056     , value_type  eta
00057     , int     after_k_iter
00058     , value_type  obj_increase_threshold
00059     , value_type  max_pos_penalty_increase
00060     , value_type  pos_to_neg_penalty_increase
00061     , value_type  incr_mult_factor )
00062   : merit_func_(merit_func), eta_(eta), after_k_iter_(after_k_iter)
00063     , obj_increase_threshold_(obj_increase_threshold)
00064     , max_pos_penalty_increase_(max_pos_penalty_increase)
00065     , pos_to_neg_penalty_increase_(pos_to_neg_penalty_increase)
00066     , incr_mult_factor_(incr_mult_factor)
00067 {}
00068 
00069 bool MeritFunc_ModifiedL1LargerSteps_AddedStep::do_step(Algorithm& _algo
00070   , poss_type step_poss, IterationPack::EDoStepType type
00071   , poss_type assoc_step_poss)
00072 {
00073   using DenseLinAlgPack::norm_inf;
00074   using DenseLinAlgPack::dot;
00075 
00076   NLPAlgo &algo = rsqp_algo(_algo);
00077   NLPAlgoState  &s    = algo.rsqp_state();
00078   
00079   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00080   std::ostream& out = algo.track().journal_out();
00081 
00082   // print step header.
00083   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00084     using IterationPack::print_algorithm_step;
00085     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00086   }
00087 
00088   MeritFuncPenaltyParams
00089     *params = dynamic_cast<MeritFuncPenaltyParams*>(&merit_func());
00090   if( !params ) {
00091     std::ostringstream omsg;
00092     omsg
00093       << "MeritFunc_ModifiedL1LargerSteps_AddedStep::do_step(...), Error "
00094       << "The class " << typeName(&merit_func()) << " does not support the "
00095       << "MeritFuncPenaltyParams iterface\n";
00096     out << omsg.str();
00097     throw std::logic_error( omsg.str() );
00098   }
00099 
00100   bool consider_modifications = s.k() >= after_k_iter();
00101   
00102   if( !consider_modifications )
00103     return true;
00104 
00105   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00106     out << "\nk = " << s.k() << " >= " << " after_k_iter = " << after_k_iter()
00107       << "\nConsidering increasing the penalty parameters ...\n";
00108   }
00109 
00110   // /////////////////////////////////////////
00111   // Get references to iteration quantities
00112 
00113   const value_type
00114     &f_k    = s.f().get_k(0),
00115     &f_kp1    = s.f().get_k(+1);
00116 
00117   const DVector
00118     &c_k    = s.c().get_k(0).cv(),
00119     &c_kp1    = s.c().get_k(+1).cv();
00120 
00121   const DVector
00122     &Gf_k   = s.Gf().get_k(0).cv(),
00123     &d_k    = s.d().get_k(0).cv();
00124 
00125   // Determining if the objective increase is sufficent.
00126 
00127   const value_type
00128     very_small  = std::numeric_limits<value_type>::min(),
00129     obj_increase = ( f_kp1 - f_k ) / ::fabs( f_kp1 + f_k + very_small );
00130   bool attempt_modifications = obj_increase >= obj_increase_threshold(); 
00131 
00132   if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00133     out << "\n(f_kp1-f_k)/|f_kp1+f_k+very_small| = " << obj_increase
00134       << ( attempt_modifications ? " >= " : " < " )
00135       << "obj_increase_threshold = " << obj_increase_threshold() << std::endl;
00136   }
00137 
00138   if( obj_increase < obj_increase_threshold() ) {
00139     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00140       out << "\nLeave the penalty parameters as they are.\n";
00141     }
00142   }
00143   else {
00144     // Compute the penalty parameters.
00145     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00146       out << "\nTry to increase the penalty parameters to allow a larger SQP step... .\n";
00147     }
00148 
00149     DVectorSlice
00150       mu = params->mu();
00151 
00152     // ///////////////////////////////////////////////////////////
00153     // Derivation of the modification to the penalty parameters
00154     //
00155     // Given the modified L1 merit function:
00156     //
00157     // phi(x) = f(x) + sum( mu(j) * |c(x)(j)|, j=1..m )
00158     // Dphi(x_k,d) = Gf_k'*d - sum( mu(j) * |c(x)(j)|, j=1..m )
00159     //
00160     // Given the armijo condition for a full step:
00161     //
00162     // phi(x_kp1) <= phi(x_k) + eta * Dphi(x_k,d)
00163     //
00164     // ->
00165     //
00166     // f_kp1 - sum(mu(j)*|c_kp1(j)|,j=1..m)
00167     //    <= f_k - sum(mu(j)*|c_k(j)|,j=1..m)
00168     //      + eta*( Gf_k'*d - sum(mu(j)*|c_k(j)|,j=1..m) )
00169     //
00170     // -> (1)
00171     //
00172     // f_kp1 - f_k - eta*Gf_k'*d <= sum( mu(j) * del(j), j=1...m )
00173     //
00174     // where:
00175     //    del(j) = (1-eta) * c_k(j) - c_kp1(j)
00176     //
00177     // Define the sets:
00178     //
00179     // DelPos = { j | del(j) >  0 }   (2.a)
00180     // DelNeg = { j | del(j) =< 0 }   (2.b)
00181     //
00182     // Define the update expresions for the penalty parameters:
00183     //
00184     // mu(j) <- mu(j) + a * ( mu_max - mu(j) ), for j <: DelPos     (3.a)
00185     //
00186     // mu(j) <- mu(j) + (a/b) * ( mu_max - mu(j) ), for j <: DelNeg (3.b)
00187     //
00188     // where:
00189     //    a < max_pos_penalty_increase ( >= 0 )
00190     //    b = pos_to_neg_penalty_increase ( >= 0 )
00191     //    mu_max = (1.0 + incr_mult_factor) * ||mu||inf
00192     //    0 < a < max_pos_penalty_increase : The length to be determined
00193     //      so that (1) can be satsified.
00194     //
00195     // The idea there is to make (1) an equality, plug (3) into it and then
00196     // solve for a.
00197     //
00198     // (1), (3) -> (4)
00199     //
00200     // a = ( f_kp1 - f_k - eta*Gf_k'*d - num_term ) / ( pos_denom_term + neg_denom_term )
00201     //
00202     // where:
00203     //    num_term = sum( mu(j) * del(j), j=1..m )
00204     //    pos_denom_term = sum( (mu_max - mu(j)) * del(j), j <: DelPos )
00205     //    neg_denom_term = (1/b) * sum( (mu_max - mu(j)) * del(j), j <: NegPos )
00206     //
00207     // If the value of a from (4) is within 0 < a < max_pos_penalty_increase
00208     // then that means that can increase the penalties
00209     // enough and satisfy (1).  If a < 0 then we would
00210     // have to decrease the penalties and we are not allowed to do this.
00211     // If a > max_pos_penalty_increase then we are not allowed to increase
00212     // the penalties enough to
00213     // satisfy (1) but it suggests that if we increase them up to (3) for
00214     // a = max_pos_penalty_increase
00215     // that we would be able to take a larger SQP step durring our linesearch.
00216 
00217     // //////////////////////////
00218     // Compute the terms in (4)
00219 
00220     const value_type
00221       mu_max = norm_inf( mu ) * (1.0 + incr_mult_factor());
00222 
00223     value_type
00224       num_term = 0.0,
00225       pos_denom_term = 0.0,
00226       neg_denom_term = 0.0;
00227 
00228     typedef std::vector<bool> del_pos_t;  // Remember the sets DelPos, DelNeg
00229     del_pos_t
00230       del_pos( mu.size() );
00231 
00232     {
00233       DVectorSlice::const_iterator
00234         mu_itr    = const_cast<const DVectorSlice&>(mu).begin();
00235       DVector::const_iterator
00236         c_k_itr   = c_k.begin(),
00237         c_kp1_itr = c_kp1.begin();
00238 
00239       del_pos_t::iterator
00240         del_pos_itr = del_pos.begin();
00241 
00242       for( ; c_k_itr != c_k.end(); ++mu_itr, ++c_k_itr, ++c_kp1_itr, ++del_pos_itr ) {
00243         TEST_FOR_EXCEPT( !(  mu_itr < const_cast<const DVectorSlice&>(mu).end()  ) );
00244         TEST_FOR_EXCEPT( !(  c_kp1_itr < c_kp1.end()  ) );
00245         TEST_FOR_EXCEPT( !(  del_pos_itr < del_pos.end()  ) );
00246 
00247         const value_type
00248           del_j = ( 1 - eta() ) * ::fabs( *c_k_itr ) - ::fabs( *c_kp1_itr );
00249           
00250         num_term += (*mu_itr) * del_j;
00251 
00252         if( del_j > 0 ) {
00253           *del_pos_itr = true;
00254           pos_denom_term += ( mu_max - (*mu_itr) ) * del_j;
00255         }
00256         else {
00257           *del_pos_itr = false;
00258           neg_denom_term += ( mu_max - (*mu_itr) ) * del_j;
00259         }
00260       }
00261       neg_denom_term /= pos_to_neg_penalty_increase();
00262     }
00263 
00264     // Compute a from (4)
00265     value_type
00266       a = ( f_kp1 - f_k - eta() * dot(Gf_k,d_k) - num_term)
00267         / ( pos_denom_term + neg_denom_term );
00268 
00269     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00270       out << "\nnum_term       = " << num_term
00271         << "\npos_denom_term = " << pos_denom_term
00272         << "\nneg_denom_term = " << neg_denom_term
00273         << "\n\na = " << a << std::endl;
00274     }
00275 
00276     if( a < 0.0 ) {
00277       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00278         out << "\na < 0 : Leave the penalty parameters alone\n";
00279       }
00280       return true;
00281     }
00282     else if( a > max_pos_penalty_increase() ) {
00283       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00284         out << "\na > max_pos_penalty_increase = " << max_pos_penalty_increase()
00285           << "\nSet a = max_pos_penalty_increase ...\n";
00286       }
00287       a = max_pos_penalty_increase();
00288     }
00289     else {
00290       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00291         out << "\n0 <= a <= max_pos_penalty_increase = " << max_pos_penalty_increase()
00292           << "\nWe should be able to take a full SQP step ...\n";
00293       }
00294     }
00295 
00296     // Update the penalty parameters using (3)
00297     {
00298       const value_type
00299         pos_step = a,
00300         neg_step = pos_step / pos_to_neg_penalty_increase(); 
00301       del_pos_t::const_iterator
00302         del_pos_itr = const_cast<const del_pos_t&>(del_pos).begin();
00303       DVectorSlice::iterator
00304         mu_itr    = mu.begin();
00305       for( ; mu_itr != mu.end(); ++del_pos_itr, ++mu_itr ) {
00306         TEST_FOR_EXCEPT( !(  del_pos_itr < const_cast<const del_pos_t&>(del_pos).end()  ) );
00307         *mu_itr = *mu_itr
00308               + (*del_pos_itr ?pos_step :neg_step) * (mu_max - (*mu_itr));
00309       }   
00310     }
00311 
00312     if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00313       out << "\nmax(|mu(j)|) = " << (*std::max_element( mu.begin(), mu.end() ))
00314         << "\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
00315           << std::endl;
00316     }
00317 
00318     if( (int)olevel >= (int)PRINT_VECTORS ) {
00319       out << "\nmu = \n" << mu;
00320     }
00321   }
00322 
00323   return true;
00324 }
00325 
00326 void MeritFunc_ModifiedL1LargerSteps_AddedStep::print_step( const Algorithm& algo
00327   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00328   , std::ostream& out, const std::string& L ) const
00329 {
00330   out
00331     << L << "*** Increase the penalty parameters for the modified L1 merit function\n"
00332     << L << "*** so as to allow for a larger SQP step.\n"
00333     << L << "default: eta                          = " << eta() << std::endl
00334     << L << "         after_k_iter                 = " << after_k_iter() << std::endl
00335     << L << "         obj_increase_threshold       = " << obj_increase_threshold() << std::endl
00336     << L << "         max_pos_penalty_increase     = " << max_pos_penalty_increase() << std::endl
00337     << L << "         pos_to_neg_penalty_increase  = " << pos_to_neg_penalty_increase() << std::endl
00338     << L << "         incr_mult_factor             = " << incr_mult_factor() << std::endl
00339     << L << "if k < after_k_iter then\n"
00340     << L << "    goto next step\n"
00341     << L << "end\n"
00342     << L << "if (f_kp1-f_k)/abs(f_kp1+f_k+very_small) >= obj_increase_threshold then\n"
00343     << L << "    mu = phi.mu()\n"
00344     << L << "    *** Try to increase to penalty parameters mu(j) to allow for a full step.\n"
00345     << L << "    mu_max = norm(mu,inf) * (1.0+incr_mult_factor)\n"
00346     << L << "    num_term = 0\n"
00347     << L << "    pos_denom_term = 0\n"
00348     << L << "    neg_denom_term = 0\n"
00349     << L << "    for j = 1 ... m\n"
00350     << L << "        del(j) = (1-eta)*abs(c_k(j)) - abs(c_kp1(k))\n"
00351     << L << "        num_term = num_term + mu(j) * del(j)\n"
00352     << L << "        if del(j) > 0 then\n"
00353     << L << "            del_pos(j) = true\n"
00354     << L << "            pos_denom_term = pos_denom_term + (mu_max - mu(j)) * del(j)\n"
00355     << L << "        else\n"
00356     << L << "            del_pos(j) = false\n"
00357     << L << "            neg_denom_term = neg_denom_term + (mu_max - mu(j)) * del(j)\n"
00358     << L << "        end\n"
00359     << L << "    end\n"
00360     << L << "    neg_denom_term = (1/pos_to_neg_penalty_increase) * neg_denom_term\n"
00361     << L << "    a = ( f_kp1 - f_k - eta * dot(Gf_k,d_k) - num_term)\n"
00362     << L << "       / ( pos_denom_term + neg_denom_term )\n"
00363     << L << "    if a < 0 then\n"
00364     << L << "        *** We can't take a larger SQP step by increasing mu(j)\n"
00365     << L << "        goto next step\n"
00366     << L << "    else if a > max_pos_penalty_increase then\n"
00367     << L << "        *** We are not allowed to increase mu(j) enough to allow a\n"
00368     << L << "        *** full SQP step but we will still increase mu(j) to take\n"
00369     << L << "        *** a hopefully larger step\n"
00370     << L << "        a = max_pos_penalty_increase\n"
00371     << L << "    else\n"
00372     << L << "        *** We can increase mu(j) and take a full SQP step\n"
00373     << L << "    end\n"
00374     << L << "    *** Increase the multipliers\n"
00375     << L << "    for j = 1...m\n"
00376     << L << "        if del_pos(j) == true then\n"
00377     << L << "            mu(j) = mu(j) + a*(mu_max - mu(j))\n"
00378     << L << "        else\n"
00379     << L << "            mu(j) = mu(j) + (a/pos_to_neg_penalty_increase)*(mu_max - mu(j))\n"
00380     << L << "        end\n"
00381     << L << "    end\n"
00382     << L << "end\n";
00383 }
00384 
00385 } // end namespace MoochoPack
00386 
00387 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines