ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <ostream>
43 #include <iomanip>
44 #include <sstream>
45 
46 #include "ConstrainedOptPack_DirectLineSearchArmQuad_Strategy.hpp"
47 #include "ConstrainedOptPack_MeritFuncCalc1D.hpp"
48 #include "check_nan_inf.h"
49 
50 namespace ConstrainedOptPack {
51 inline value_type min(value_type v1, value_type v2) {
52  return (v1 < v2) ? v1 : v2;
53 }
54 inline value_type max(value_type v1, value_type v2) {
55  return (v1 > v2) ? v1 : v2;
56 }
57 }
58 
60  int max_iter
61  ,value_type eta
62  ,value_type min_frac
63  ,value_type max_frac
64  ,bool max_out_iter
65  )
66  :max_iter_(max_iter)
67  ,eta_(eta)
68  ,min_frac_(min_frac)
69  ,max_frac_(max_frac)
70  ,max_out_iter_(max_out_iter)
71 {}
72 
74 {
75  max_iter_ = max_iter;
76 }
77 
79 {
80  return max_iter_;
81 }
82 
84 {
85  return num_iter_;
86 }
87 
89  const MeritFuncCalc1D& phi, value_type phi_k
90  , value_type* alpha_k, value_type* phi_kp1
91  , std::ostream* out )
92 {
93  using std::setw;
94  using std::endl;
95 
96  if(*alpha_k < 0.0) {
97  throw std::logic_error( "DirectLineSearchArmQuad_Strategy::do_line_search(): "
98  "alpha_k can't start out less than 0.0" );
99  }
100 
101  validate_parameters();
102 
103  int w = 20;
104  int prec = 8;
105  int prec_saved;
106  if(out) {
107  prec_saved = out->precision();
108  *out << std::setprecision(prec)
109  << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
110  }
111 
112  // Loop initialization (technically the first iteration)
113 
114  const value_type Dphi_k = phi.deriv();
115 
116  // output header
117  if(out) {
118  if(max_out_iter())
119  *out << "\nmax_out_iter == true, maxing out the number of iterations!";
120  *out << "\nDphi_k = " << Dphi_k
121  << "\nphi_k = " << phi_k << "\n\n"
122  << setw(5) << "itr"
123  << setw(w) << "alpha_k"
124  << setw(w) << "phi_kp1"
125  << setw(w) << "phi_kp1-frac_phi\n"
126  << setw(5) << "----"
127  << setw(w) << "----------------"
128  << setw(w) << "----------------"
129  << setw(w) << "----------------\n";
130  }
131 
132  // Check that this is a decent direction
134  "DirectLineSearchArmQuad_Strategy::do_line_search(): "
135  "The given d_k is not a descent direction for the given "
136  "phi (phi.deriv() is >= 0)" );
137 
138  // keep memory of the best value
139  value_type best_alpha = *alpha_k, best_phi = *phi_kp1;
140 
141  // Perform linesearch.
142  bool success = false;
143  for( num_iter_ = 0; num_iter_ < max_iter(); ++num_iter_ ) {
144 
145  // Print out this iteration.
146 
147  value_type frac_phi = phi_k + eta() * (*alpha_k) * Dphi_k;
148  if(out)
149  *out << setw(5) << num_iter_
150  << setw(w) << *alpha_k
151  << setw(w) << *phi_kp1
152  << setw(w) << ((*phi_kp1)-frac_phi) << endl;
153 
154  // Check that this is a valid number.
155  if( RTOp_is_nan_inf( *phi_kp1 ) ) {
156  // Cut back the step to min_frac * alpha_k
157  *alpha_k = min_frac()*(*alpha_k);
158  best_alpha = 0.0;
159  best_phi = phi_k;
160  }
161  else {
162 
163  // Armijo condition
164  if( *phi_kp1 < frac_phi ) {
165  // We have found an acceptable point
166  success = true;
167  if( !max_out_iter() || ( max_out_iter() && num_iter_ == max_iter() - 1 ) )
168  break; // get out of the loop
169  }
170 
171  // Select a new alpha to try:
172  // alpha_k = ( min_frac*alpha_k <= quadratic interpolation <= max_frac*alpha_k )
173 
174  // Quadratic interpolation of alpha_k that minimizes phi.
175  // We know the values of phi at the initail point and alpha_k and
176  // the derivate of phi w.r.t alpha at the initial point and
177  // that's enough information for a quandratic interpolation.
178 
179  value_type alpha_quad = ( -0.5 * Dphi_k * (*alpha_k) * (*alpha_k) )
180  / ( (*phi_kp1) - phi_k - (*alpha_k) * Dphi_k );
181 
182  *alpha_k = min( max(min_frac()*(*alpha_k),alpha_quad), max_frac()*(*alpha_k) );
183 
184  }
185 
186  // Evaluate the point
187 
188  *phi_kp1 = phi(*alpha_k);
189 
190  // Save the best point found
191  if(*phi_kp1 < best_phi) {
192  best_phi = *phi_kp1;
193  best_alpha = *alpha_k;
194  }
195 
196  }
197 
198  // Be nice and reset the precision
199  if(out) {
200  out->precision(prec_saved);
201  }
202 
203  if( success ) {
204  return true;
205  }
206 
207  // Line search failure. Return the best point found and let the
208  // client decide what to do.
209  *alpha_k = best_alpha;
210  *phi_kp1 = phi(best_alpha); // Make this the last call to phi(x)
211  return false;
212 }
213 
215  std::ostream& out, const std::string& L) const
216 {
217  out
218  << L << "*** start line search using the Armijo cord test and quadratic interpolation of alpha\n"
219  << L << "default: max_ls_iter = " << max_iter() << std::endl
220  << L << " eta = " << eta() << std::endl
221  << L << " min_frac = " << min_frac() << std::endl
222  << L << " max_frac = " << max_frac() << std::endl
223  << L << " max_out_iter = " << max_out_iter() << std::endl
224  << L << "Dphi_k = phi.deriv()\n"
225  << L << "if Dphi_k >= 0\n"
226  << L << " throw not_descent_direction\n"
227  << L << " end line search\n"
228  << L << "end\n"
229  << L << "best_alpha = alpha_k\n"
230  << L << "best_phi = phi_kp1\n"
231  << L << "for num_iter = 0... max_ls_iter\n"
232  << L << " frac_phi = phi_k + eta * alpha_k * Dphi_k\n"
233  << L << " print iteration\n"
234  << L << " if phi_kp1 is not a valid number then\n"
235  << L << " *** Cut back the step so the NLP's functions\n"
236  << L << " *** will hopefully be defined.\n"
237  << L << " alpha_k = min_frac * alpha_k\n"
238  << L << " best_alpha = 0\n"
239  << L << " best_phi = phi_k\n"
240  << L << " else\n"
241  << L << " if ( phi_kp1 < frac_phi ) then\n"
242  << L << " *** We have found an acceptable point\n"
243  << L << " if( !max_out_iter || max_out_iter && num_iter == max_ls_iter - 1 ) )\n"
244  << L << " end line search\n"
245  << L << " end\n"
246  << L << " end\n"
247  << L << " *** Use a quadratic interpoation to minimize phi(alpha)\n"
248  << L << " alpha_quad = (-0.5 * Dphi_k * alpha_k^2) / ( phi_kp1 - phi_k - alpha_k*Dphi_k )\n"
249  << L << " alpha_k = min( max( min_frac*alpha_k, alpha_quad ), max_frac*alpha_k )\n"
250  << L << " end\n"
251  << L << " phi_kp1 = phi(alpha_k)\n"
252  << L << " if phi_kp1 < best_phi\n"
253  << L << " best_phi = phi_kp1\n"
254  << L << " best_alpha = alpha_k\n"
255  << L << " end\n"
256  << L << "end\n"
257  << L << "*** If you get there the line search failed.\n"
258  << L << "alpha_k = best_alpha\n"
259  << L << "phi_kp1 = phi(alpha_k)\n"
260  << L << "linesearch_failure = true\n";
261 }
262 
263 void ConstrainedOptPack::DirectLineSearchArmQuad_Strategy::validate_parameters() const
264 {
265  if( eta() < 0.0 || 1.0 < eta() ) {
266  std::ostringstream omsg;
267  omsg
268  << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
269  << "Error, eta = " << eta() << " is not in the range [0, 1]";
270  throw std::invalid_argument( omsg.str() );
271  }
272  if( min_frac() < 0.0 || 1.0 < min_frac() ) {
273  std::ostringstream omsg;
274  omsg
275  << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
276  << "Error, min_frac = " << min_frac() << " is not in the range [0, 1]";
277  throw std::invalid_argument( omsg.str() );
278  }
279  if( max_frac() < 0.0 || ( !max_out_iter() && 1.0 < max_frac() ) ) {
280  std::ostringstream omsg;
281  omsg
282  << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
283  << "Error, max_frac = " << max_frac() << " is not in the range [0, 1]";
284  throw std::invalid_argument( omsg.str() );
285  }
286  if( max_frac() < min_frac() ) {
287  std::ostringstream omsg;
288  omsg
289  << "DirectLineSearchArmQuad_Strategy::validate_parameters() : "
290  << "Error, max_frac = " << max_frac()
291  << " < min_frac = " << min_frac();;
292  throw std::invalid_argument( omsg.str() );
293  }
294 }
Abstracts a 1D merit function {abstract}.
bool do_line_search(const MeritFuncCalc1D &phi, value_type phi_k, value_type *alpha_k, value_type *phi_kp1, std::ostream *out)
Performs the following line search:
void print_algorithm(std::ostream &out, const std::string &leading_str) const
Thrown if the direction vector d_k is not a descent direction for the merit funciton.
virtual value_type deriv() const =0
Return the derivative of the merit function at alpha = 0.
DirectLineSearchArmQuad_Strategy(int max_iter=20, value_type eta=1.0e-4, value_type min_frac=0.1, value_type max_frac=0.5, bool max_out_iter=false)
Constructs with default settings.