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

Generated on Thu Sep 18 12:35:17 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1