MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_MeritFunc_ModifiedL1LargerSteps_AddedStep.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <ostream>
45 #include <typeinfo>
46 
47 #include "MoochoPack_MeritFunc_ModifiedL1LargerSteps_AddedStep.hpp"
48 #include "MoochoPack_moocho_algo_conversion.hpp"
49 #include "IterationPack_print_algorithm_step.hpp"
50 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
51 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
52 #include "ConstrainedOptPack/src/VectorWithNorms.h"
53 #include "DenseLinAlgPack_DVectorOp.hpp"
54 #include "DenseLinAlgPack_DVectorClass.hpp"
55 #include "DenseLinAlgPack_DVectorOut.hpp"
56 
57 namespace {
58 
59 typedef MoochoPack::value_type value_type;
60 inline value_type max(value_type v1, value_type v2)
61 { return (v1 > v2) ? v1 : v2; }
62 
63 }
64 
65 namespace MoochoPack {
66 
67 MeritFunc_ModifiedL1LargerSteps_AddedStep::MeritFunc_ModifiedL1LargerSteps_AddedStep(
68  const merit_func_ptr_t& merit_func
69  , value_type eta
70  , int after_k_iter
71  , value_type obj_increase_threshold
72  , value_type max_pos_penalty_increase
73  , value_type pos_to_neg_penalty_increase
74  , value_type incr_mult_factor )
75  : merit_func_(merit_func), eta_(eta), after_k_iter_(after_k_iter)
76  , obj_increase_threshold_(obj_increase_threshold)
77  , max_pos_penalty_increase_(max_pos_penalty_increase)
78  , pos_to_neg_penalty_increase_(pos_to_neg_penalty_increase)
79  , incr_mult_factor_(incr_mult_factor)
80 {}
81 
83  , poss_type step_poss, IterationPack::EDoStepType type
84  , poss_type assoc_step_poss)
85 {
86  using DenseLinAlgPack::norm_inf;
88 
89  NLPAlgo &algo = rsqp_algo(_algo);
90  NLPAlgoState &s = algo.rsqp_state();
91 
92  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
93  std::ostream& out = algo.track().journal_out();
94 
95  // print step header.
96  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
97  using IterationPack::print_algorithm_step;
98  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
99  }
100 
101  MeritFuncPenaltyParams
102  *params = dynamic_cast<MeritFuncPenaltyParams*>(&merit_func());
103  if( !params ) {
104  std::ostringstream omsg;
105  omsg
106  << "MeritFunc_ModifiedL1LargerSteps_AddedStep::do_step(...), Error "
107  << "The class " << typeName(&merit_func()) << " does not support the "
108  << "MeritFuncPenaltyParams iterface\n";
109  out << omsg.str();
110  throw std::logic_error( omsg.str() );
111  }
112 
113  bool consider_modifications = s.k() >= after_k_iter();
114 
115  if( !consider_modifications )
116  return true;
117 
118  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
119  out << "\nk = " << s.k() << " >= " << " after_k_iter = " << after_k_iter()
120  << "\nConsidering increasing the penalty parameters ...\n";
121  }
122 
123  // /////////////////////////////////////////
124  // Get references to iteration quantities
125 
126  const value_type
127  &f_k = s.f().get_k(0),
128  &f_kp1 = s.f().get_k(+1);
129 
130  const DVector
131  &c_k = s.c().get_k(0).cv(),
132  &c_kp1 = s.c().get_k(+1).cv();
133 
134  const DVector
135  &Gf_k = s.Gf().get_k(0).cv(),
136  &d_k = s.d().get_k(0).cv();
137 
138  // Determining if the objective increase is sufficent.
139 
140  const value_type
141  very_small = std::numeric_limits<value_type>::min(),
142  obj_increase = ( f_kp1 - f_k ) / ::fabs( f_kp1 + f_k + very_small );
143  bool attempt_modifications = obj_increase >= obj_increase_threshold();
144 
145  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
146  out << "\n(f_kp1-f_k)/|f_kp1+f_k+very_small| = " << obj_increase
147  << ( attempt_modifications ? " >= " : " < " )
148  << "obj_increase_threshold = " << obj_increase_threshold() << std::endl;
149  }
150 
151  if( obj_increase < obj_increase_threshold() ) {
152  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
153  out << "\nLeave the penalty parameters as they are.\n";
154  }
155  }
156  else {
157  // Compute the penalty parameters.
158  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
159  out << "\nTry to increase the penalty parameters to allow a larger SQP step... .\n";
160  }
161 
162  DVectorSlice
163  mu = params->mu();
164 
165  // ///////////////////////////////////////////////////////////
166  // Derivation of the modification to the penalty parameters
167  //
168  // Given the modified L1 merit function:
169  //
170  // phi(x) = f(x) + sum( mu(j) * |c(x)(j)|, j=1..m )
171  // Dphi(x_k,d) = Gf_k'*d - sum( mu(j) * |c(x)(j)|, j=1..m )
172  //
173  // Given the armijo condition for a full step:
174  //
175  // phi(x_kp1) <= phi(x_k) + eta * Dphi(x_k,d)
176  //
177  // ->
178  //
179  // f_kp1 - sum(mu(j)*|c_kp1(j)|,j=1..m)
180  // <= f_k - sum(mu(j)*|c_k(j)|,j=1..m)
181  // + eta*( Gf_k'*d - sum(mu(j)*|c_k(j)|,j=1..m) )
182  //
183  // -> (1)
184  //
185  // f_kp1 - f_k - eta*Gf_k'*d <= sum( mu(j) * del(j), j=1...m )
186  //
187  // where:
188  // del(j) = (1-eta) * c_k(j) - c_kp1(j)
189  //
190  // Define the sets:
191  //
192  // DelPos = { j | del(j) > 0 } (2.a)
193  // DelNeg = { j | del(j) =< 0 } (2.b)
194  //
195  // Define the update expresions for the penalty parameters:
196  //
197  // mu(j) <- mu(j) + a * ( mu_max - mu(j) ), for j <: DelPos (3.a)
198  //
199  // mu(j) <- mu(j) + (a/b) * ( mu_max - mu(j) ), for j <: DelNeg (3.b)
200  //
201  // where:
202  // a < max_pos_penalty_increase ( >= 0 )
203  // b = pos_to_neg_penalty_increase ( >= 0 )
204  // mu_max = (1.0 + incr_mult_factor) * ||mu||inf
205  // 0 < a < max_pos_penalty_increase : The length to be determined
206  // so that (1) can be satsified.
207  //
208  // The idea there is to make (1) an equality, plug (3) into it and then
209  // solve for a.
210  //
211  // (1), (3) -> (4)
212  //
213  // a = ( f_kp1 - f_k - eta*Gf_k'*d - num_term ) / ( pos_denom_term + neg_denom_term )
214  //
215  // where:
216  // num_term = sum( mu(j) * del(j), j=1..m )
217  // pos_denom_term = sum( (mu_max - mu(j)) * del(j), j <: DelPos )
218  // neg_denom_term = (1/b) * sum( (mu_max - mu(j)) * del(j), j <: NegPos )
219  //
220  // If the value of a from (4) is within 0 < a < max_pos_penalty_increase
221  // then that means that can increase the penalties
222  // enough and satisfy (1). If a < 0 then we would
223  // have to decrease the penalties and we are not allowed to do this.
224  // If a > max_pos_penalty_increase then we are not allowed to increase
225  // the penalties enough to
226  // satisfy (1) but it suggests that if we increase them up to (3) for
227  // a = max_pos_penalty_increase
228  // that we would be able to take a larger SQP step durring our linesearch.
229 
230  // //////////////////////////
231  // Compute the terms in (4)
232 
233  const value_type
234  mu_max = norm_inf( mu ) * (1.0 + incr_mult_factor());
235 
236  value_type
237  num_term = 0.0,
238  pos_denom_term = 0.0,
239  neg_denom_term = 0.0;
240 
241  typedef std::vector<bool> del_pos_t; // Remember the sets DelPos, DelNeg
242  del_pos_t
243  del_pos( mu.size() );
244 
245  {
246  DVectorSlice::const_iterator
247  mu_itr = const_cast<const DVectorSlice&>(mu).begin();
248  DVector::const_iterator
249  c_k_itr = c_k.begin(),
250  c_kp1_itr = c_kp1.begin();
251 
252  del_pos_t::iterator
253  del_pos_itr = del_pos.begin();
254 
255  for( ; c_k_itr != c_k.end(); ++mu_itr, ++c_k_itr, ++c_kp1_itr, ++del_pos_itr ) {
256  TEUCHOS_TEST_FOR_EXCEPT( !( mu_itr < const_cast<const DVectorSlice&>(mu).end() ) );
257  TEUCHOS_TEST_FOR_EXCEPT( !( c_kp1_itr < c_kp1.end() ) );
258  TEUCHOS_TEST_FOR_EXCEPT( !( del_pos_itr < del_pos.end() ) );
259 
260  const value_type
261  del_j = ( 1 - eta() ) * ::fabs( *c_k_itr ) - ::fabs( *c_kp1_itr );
262 
263  num_term += (*mu_itr) * del_j;
264 
265  if( del_j > 0 ) {
266  *del_pos_itr = true;
267  pos_denom_term += ( mu_max - (*mu_itr) ) * del_j;
268  }
269  else {
270  *del_pos_itr = false;
271  neg_denom_term += ( mu_max - (*mu_itr) ) * del_j;
272  }
273  }
274  neg_denom_term /= pos_to_neg_penalty_increase();
275  }
276 
277  // Compute a from (4)
278  value_type
279  a = ( f_kp1 - f_k - eta() * dot(Gf_k,d_k) - num_term)
280  / ( pos_denom_term + neg_denom_term );
281 
282  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
283  out << "\nnum_term = " << num_term
284  << "\npos_denom_term = " << pos_denom_term
285  << "\nneg_denom_term = " << neg_denom_term
286  << "\n\na = " << a << std::endl;
287  }
288 
289  if( a < 0.0 ) {
290  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
291  out << "\na < 0 : Leave the penalty parameters alone\n";
292  }
293  return true;
294  }
295  else if( a > max_pos_penalty_increase() ) {
296  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
297  out << "\na > max_pos_penalty_increase = " << max_pos_penalty_increase()
298  << "\nSet a = max_pos_penalty_increase ...\n";
299  }
300  a = max_pos_penalty_increase();
301  }
302  else {
303  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
304  out << "\n0 <= a <= max_pos_penalty_increase = " << max_pos_penalty_increase()
305  << "\nWe should be able to take a full SQP step ...\n";
306  }
307  }
308 
309  // Update the penalty parameters using (3)
310  {
311  const value_type
312  pos_step = a,
313  neg_step = pos_step / pos_to_neg_penalty_increase();
314  del_pos_t::const_iterator
315  del_pos_itr = const_cast<const del_pos_t&>(del_pos).begin();
316  DVectorSlice::iterator
317  mu_itr = mu.begin();
318  for( ; mu_itr != mu.end(); ++del_pos_itr, ++mu_itr ) {
319  TEUCHOS_TEST_FOR_EXCEPT( !( del_pos_itr < const_cast<const del_pos_t&>(del_pos).end() ) );
320  *mu_itr = *mu_itr
321  + (*del_pos_itr ?pos_step :neg_step) * (mu_max - (*mu_itr));
322  }
323  }
324 
325  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
326  out << "\nmax(|mu(j)|) = " << (*std::max_element( mu.begin(), mu.end() ))
327  << "\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
328  << std::endl;
329  }
330 
331  if( (int)olevel >= (int)PRINT_VECTORS ) {
332  out << "\nmu = \n" << mu;
333  }
334  }
335 
336  return true;
337 }
338 
340  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
341  , std::ostream& out, const std::string& L ) const
342 {
343  out
344  << L << "*** Increase the penalty parameters for the modified L1 merit function\n"
345  << L << "*** so as to allow for a larger SQP step.\n"
346  << L << "default: eta = " << eta() << std::endl
347  << L << " after_k_iter = " << after_k_iter() << std::endl
348  << L << " obj_increase_threshold = " << obj_increase_threshold() << std::endl
349  << L << " max_pos_penalty_increase = " << max_pos_penalty_increase() << std::endl
350  << L << " pos_to_neg_penalty_increase = " << pos_to_neg_penalty_increase() << std::endl
351  << L << " incr_mult_factor = " << incr_mult_factor() << std::endl
352  << L << "if k < after_k_iter then\n"
353  << L << " goto next step\n"
354  << L << "end\n"
355  << L << "if (f_kp1-f_k)/abs(f_kp1+f_k+very_small) >= obj_increase_threshold then\n"
356  << L << " mu = phi.mu()\n"
357  << L << " *** Try to increase to penalty parameters mu(j) to allow for a full step.\n"
358  << L << " mu_max = norm(mu,inf) * (1.0+incr_mult_factor)\n"
359  << L << " num_term = 0\n"
360  << L << " pos_denom_term = 0\n"
361  << L << " neg_denom_term = 0\n"
362  << L << " for j = 1 ... m\n"
363  << L << " del(j) = (1-eta)*abs(c_k(j)) - abs(c_kp1(k))\n"
364  << L << " num_term = num_term + mu(j) * del(j)\n"
365  << L << " if del(j) > 0 then\n"
366  << L << " del_pos(j) = true\n"
367  << L << " pos_denom_term = pos_denom_term + (mu_max - mu(j)) * del(j)\n"
368  << L << " else\n"
369  << L << " del_pos(j) = false\n"
370  << L << " neg_denom_term = neg_denom_term + (mu_max - mu(j)) * del(j)\n"
371  << L << " end\n"
372  << L << " end\n"
373  << L << " neg_denom_term = (1/pos_to_neg_penalty_increase) * neg_denom_term\n"
374  << L << " a = ( f_kp1 - f_k - eta * dot(Gf_k,d_k) - num_term)\n"
375  << L << " / ( pos_denom_term + neg_denom_term )\n"
376  << L << " if a < 0 then\n"
377  << L << " *** We can't take a larger SQP step by increasing mu(j)\n"
378  << L << " goto next step\n"
379  << L << " else if a > max_pos_penalty_increase then\n"
380  << L << " *** We are not allowed to increase mu(j) enough to allow a\n"
381  << L << " *** full SQP step but we will still increase mu(j) to take\n"
382  << L << " *** a hopefully larger step\n"
383  << L << " a = max_pos_penalty_increase\n"
384  << L << " else\n"
385  << L << " *** We can increase mu(j) and take a full SQP step\n"
386  << L << " end\n"
387  << L << " *** Increase the multipliers\n"
388  << L << " for j = 1...m\n"
389  << L << " if del_pos(j) == true then\n"
390  << L << " mu(j) = mu(j) + a*(mu_max - mu(j))\n"
391  << L << " else\n"
392  << L << " mu(j) = mu(j) + (a/pos_to_neg_penalty_increase)*(mu_max - mu(j))\n"
393  << L << " end\n"
394  << L << " end\n"
395  << L << "end\n";
396 }
397 
398 } // end namespace MoochoPack
399 
400 #endif // 0
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
value_type max_element(const Vector &v)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string typeName(const T &t)