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