Interior-point-optimisation  1.0-1
Interior-pointoptimisationlibrary
NewtonDescent.cc
Go to the documentation of this file.
1 /*
2  * $Id: NewtonDescent.cc 215 2014-11-04 08:58:35Z jdl3 $
3  * Copyright (C) 2011, 2012, 2013 John D Lamb
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or (at
8  * your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful, but
11  * WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18  */
19 
20 #ifdef HAVE_CONFIG_H
21 # include<config.h>
22 #endif
23 
24 #include"NewtonDescent.hpp"
25 
26 using namespace ipo::detail;
27 
28 namespace {
29  double constexpr gradTolDefault
30  { std::exp( std::log( std::numeric_limits<double>::epsilon() ) / 3 ) };
31  double constexpr stepTolDefault { gradTolDefault * gradTolDefault };
32 }
33 
36 
39 
41  setOutputStream( nullptr );
42  setPhaseI( false );
43  setGradTol( gradTolDefault );
44  setStepTol( stepTolDefault );
45  setIterMax( 10000 );
46  setf_typ( 1 );
47  setx_typ( gsl::vector() );
48  setVectorDigits( 0 );
49 }
50 
53  gsl::matrix const equalityConstraintsMatrix,
54  gsl::vector const equalityConstraintsBounds )
55  : equalityConstraintsMatrix { equalityConstraintsMatrix },
56  function( function ), de( de ), lineSearch { function }
57 {
58  // set function scale
59  functionScale = 1;
60 
61  // Check equalityConstraintsObjects
62  if( equalityConstraintsMatrix ){
63  size_t const SIZE1 { equalityConstraintsMatrix.size1() };
64  size_t const SIZE2 { equalityConstraintsMatrix.size2() };
65  if( SIZE1 > 0 and SIZE2 > 0 ){
66  // non-null matrix
67  if( (not equalityConstraintsBounds) or 0 == equalityConstraintsBounds.size() )
68  throw IPOE( "ipo::detail::NewtonDescent::NewtonDescent(): "
69  "null equality bounds vector with non-null equality constraints matrix." );
70  else if( equalityConstraintsBounds.size() != SIZE1 )
71  throw IPOE( "ipo::detail::NewtonDescent::NewtonDescent(): "
72  "equality bounds vector equality constraints matrix sizes do not match." );
73  // Allocate matrices for singular value decomposition
74  A = gsl::matrix( SIZE1 + SIZE2, SIZE1 + SIZE2 );
75  V = gsl::matrix( SIZE1 + SIZE2, SIZE1 + SIZE2 );
76  S = gsl::vector( SIZE1 + SIZE2 );
77  work = gsl::vector( SIZE1 + SIZE2 );
78  b = gsl::vector( SIZE1 + SIZE2 );
79  for( size_t i { 0 }; i < SIZE1; ++i ) b[i + SIZE2] = equalityConstraintsBounds[i];
80  } else {
81  // Check vector is also null
82  if( equalityConstraintsBounds and equalityConstraintsBounds.size() > 0 )
83  throw IPOE( "ipo::detail::NewtonDescent::NewtonDescent(): "
84  "non-null equality bounds vector with null equality constraints matrix." );
85  }
86  } else {
87  // Check vector is also null
88  if( equalityConstraintsBounds and equalityConstraintsBounds.size() > 0 )
89  throw IPOE( "ipo::detail::NewtonDescent::NewtonDescent(): "
90  "non-null equality bounds vector with null equality constraints matrix." );
91  }
92 }
93 
94 double
95 NewtonDescent::operator()( gsl::vector& x, double objectiveValue ){
96  /* initialise various values */
97  size_t const SIZE { x.size() };
98  gsl::vector delta( SIZE );
99  gsl::vector last_x; /* deliberately set initially to evaluate to false
100  * so that stopping condition does not attempt to use it. */
101 
103  auto const outputStream = parameters.getOutputStream();
104  Stopping stopping( *this ); /* Stopping conditions */
105 
106  size_t iteration { 0 };
107  objectiveValue = function( x );
108  size_t const& digits { parameters.getVectorDigits() };
109  if( outputStream != nullptr ){
110  *outputStream << "Newton descent:" << std::endl;
111  *outputStream << "step " << ( digits > 0 ? "x " : "") << "f(x)" << std::endl;
112  *outputStream << " " << std::setw( 4 ) << iteration++ << " ";
113  if( digits > 0 ){
114  *outputStream << "( ";
115  for( size_t i = 0; i < SIZE; ++i ){
116  *outputStream << std::setw( 4 + digits ) << std::setprecision( digits ) << std::fixed
117  << x.get( i ) << " ";
118  }
119  *outputStream << ")";
120  }
121  size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
122  *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
123  << objectiveValue * functionScale << std::endl;
124  }
125 
126  /* Record final entry of vector for phase 1 */
127  double lastEntry { *x.rbegin() };
128  /* Main loop */
129  for( de.setVector( x ); stopping( iteration, objectiveValue, de.gradient(), x, last_x );
130  de.setVector( x ) ){
131  last_x = x.clone(); // Now copy x.
132 
133  de.gradient().scale( -1 ); /* negate gradient */
134 
135  gsl::matrix hessian { de.hessian() };
137  and equalityConstraintsMatrix.size2() > 0 ){
138  // We use a singular value decomposition
139  size_t const SIZE1 { equalityConstraintsMatrix.size1() };
140  size_t const SIZE2 { equalityConstraintsMatrix.size2() };
141  // Set up matrix A
142  for( size_t i { 0 }; i < SIZE2; ++i ){
143  for( size_t j { i }; j < SIZE2; ++j ){
144  double const d { hessian.get( i, j ) };
145  A.set( i, j, d );
146  if( j > i ) A.set( j, i, d );
147  }
148  for( size_t j { 0 }; j < SIZE1; ++j ){
149  double const d { equalityConstraintsMatrix.get( j, i ) };
150  A.set( SIZE2 + j, i, d );
151  A.set( i, SIZE2 + j, d );
152  }
153  }
154  for( size_t i { SIZE2 }; i < SIZE1 + SIZE2; ++i )
155  for( size_t j { i }; j < SIZE1 + SIZE2; ++j ){
156  A.set( i, j, 0 );
157  if( j > i ) A.set( j, i, 0 );
158  }
159  // Find singular value decomposition
160  auto handler = gsl::exception::set_handler_gsl_exceptions();
161  try {
162  gsl::linalg::SV_decomp( A, V, S, work );
163  // Use to solve to find delta
164  auto const& g = de.gradient();
165  for( size_t i { 0 }; i < SIZE2; ++i ) b[i] = g[i]; // rest of b is fixed as equality
166  for( size_t i { SIZE2 }; i < SIZE1 + SIZE2; ++i ) b[i] = 0;
167  gsl::linalg::SV_solve( A, V, S, b, work ); // constraints bounds in constructor
168  // copy to delta
169  for( size_t i { 0 }; i < SIZE2; ++i ) delta[i] = work[i];
170  } catch( gsl::exception& e ){
171  // FIXME: set a last stopping condition
172  return 0; // Singular value decomposition failed.
173  } catch( ... ) {
174  gsl::exception::set_handler( handler );
175  throw IPOE( "ipo::detail::NewtonDescent::operator(): "
176  "unknown exception encountered." );
177  }
178  gsl::exception::set_handler( handler );
179  } else {
180  // We can use a modified Cholesky decomposition
181  auto choleDecomp = ipo::detail::modelHessian( hessian );
182  gsl::linalg::cholesky_solve( choleDecomp, de.gradient(), delta );
183  }
184 
185  /* restore gradient */
186  de.gradient().scale( -1 ); /* negate gradient */
187  /* line search */
188  double const ls { lineSearch( x, delta, de.value(), de.gradient() ) };
189  objectiveValue = function( x );
190 
191  // Possibly put some output
192  if( outputStream != nullptr ){
193  *outputStream << " " << std::setw( 4 ) << iteration << " ";
194  if( digits > 0 ){
195  *outputStream << "( ";
196  for( size_t i { 0 }; i < SIZE; ++i ){
197  *outputStream << std::setw( 4 + digits ) << std::setprecision( digits ) << std::fixed
198  << x.get( i ) << " ";
199  }
200  *outputStream << ")";
201  }
202  size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
203  *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
204  << objectiveValue * functionScale << std::endl;
205  }
206  // Early exit possible if phase1 is true
207  double const xLastEntry { *x.rbegin() };
208  if( parameters.getPhaseI() and xLastEntry < 0 ) break;
209  // Break loop if no linesearch improvement
210  if( ls < std::numeric_limits<double>::epsilon() ) break;
211  // Break loop in phase1 if no improvement in last entry
212  if( parameters.getPhaseI() ){
213  if( xLastEntry < lastEntry ) lastEntry = xLastEntry; // force improvement
214  else break; // or force t to increase
215  }
216  ++iteration; // Finally, update iteration
217  } // end for
218 
219  // Possibly put some output
220  if( outputStream != nullptr )
221  *outputStream << "Newton descent finished." << std::endl << std::endl;
222  return objectiveValue;
223 }
224 
225 void
226 NewtonDescent::Parameters::setOutputStream( std::ostream* outputStream ){
227  this->outputStream = outputStream;
228 }
229 
230 std::ostream*
232  return outputStream;
233 }
234 
236  : newtonDescent( newtonDescent ){}
237 
238 
239 bool
240 NewtonDescent::Stopping::operator()( size_t const iteration, double const fValue,
241  gsl::vector const& gradient, gsl::vector const& vector,
242  gsl::vector const& lastVector ) const {
243  // Aliases
244  gsl::vector const& x_typ
245  { const_cast<NewtonDescent&>( newtonDescent ).getParameters().getx_typ() };
246  double const& f_typ
247  { const_cast<NewtonDescent&>( newtonDescent ).getParameters().getf_typ() };
248 
249  // Test stopping conditions one by one
250 
251  // First, test number of iterations
252  if( iteration > const_cast<NewtonDescent&>( newtonDescent ).getParameters().getIterMax() ){
253  newtonDescent.lastStoppingCondition = NewtonDescent::StoppingCondition::iterMax;
254  return false;
255  }
256 
257  // Second, test gradient
258  if( [&]() -> bool { // Return false if condition is false; else continue
259  double const GRADTOL
260  { const_cast<NewtonDescent&>( newtonDescent ).getParameters().getGradTol()
261  * std::max( f_typ, fValue ) };
262  for( size_t i { 0 }; i < vector.size(); ++i ){
263  double const n { std::fabs( gradient.get( i ) )
264  * std::max( std::fabs( vector.get( i ) ), x_typ ?
265  std::fabs( x_typ.get( i ) ) : 1 ) };
266  if( n > GRADTOL ) return true;
267  }
268  return false;
269  }() == false ){
270  newtonDescent.lastStoppingCondition = NewtonDescent::StoppingCondition::gradTol;
271  return false;
272  }
273 
274  // Third, test relative change in vector (if lastVector evaluates to true)
275  double N { 0 };
276  double D { 0 };
277  if( lastVector ){ // Typically evaluates to false at first iteration
278  for( size_t i { 0 }; i < vector.size(); ++i ){
279  double const xi { vector.get( i ) };
280  double const n { std::fabs( xi - lastVector.get( i ) ) };
281  double const d
282  { std::max( std::fabs( xi ), std::fabs( x_typ ? x_typ.get( i ) : 1 ) ) };
283  N = std::max( n, N );
284  D = std::max( d, D );
285  if( n > d * const_cast<NewtonDescent&>( newtonDescent ).getParameters().getStepTol() )
286  return true;
287  }
288  newtonDescent.lastStoppingCondition = NewtonDescent::StoppingCondition::stepTol;
289  return false;
290  }
291  return true;
292 }
293 
296  return lastStoppingCondition;
297 }
void setStepTol(double const stepTol)
Set the value of stepTol.
void setx_typ(gsl::vector x_typ)
Set the value of typical value of argument.
double functionScale
Scaling factor for output of function values.
double operator()(gsl::vector &x, double objectiveValue)
The function object operator.
double getf_typ() const
Get the value of typical value of function.
void setOutputStream(std::ostream *outputStream)
Set outputStream.
Newton descent method.
Parameters()
Constructor: sets default values.
NewtonDescent(ipo_function::Function &function, ipo_function::DerivativesEstimates &de, gsl::matrix const equalityConstraintsMatrix=gsl::matrix(), gsl::vector const equalityConstraintsBounds=gsl::vector())
The constructor needs a Function and a DerivativesEstimates (which may be identical).
Stopping(NewtonDescent const &newtonDescent)
Constructor.
gsl::vector getx_typ() const
Get the value of typical value of argument.
gsl::vector work
This vector is used in singular value decompostion and its size is set in constructor if it is needed...
void setVectorDigits(size_t const vectorDigits)
Set the number of significant digits of entries of vector to show in each Newton step.
bool operator()(size_t const iteration, double const fValue, gsl::vector const &gradient, gsl::vector const &vector, gsl::vector const &lastVector) const
Test the stopping conditions.
virtual gsl::matrix hessian() const
virtual gsl::vector gradient() const
gsl::matrix const equalityConstraintsMatrix
This must be set in the constructor.
StoppingCondition getLastStoppingCondition() const
Get Last stopping condition.
virtual double value() const
Base class for derivative estimates with Hessian.
size_t getVectorDigits() const
Get the value of vectorDigits.
void setIterMax(size_t const iterMax)
Set the value of iterMax.
Parameters & getParameters()
Get parameters by reference.
Definition: LineSearch.cc:39
void setFunctionScale(double const functionScale)
Set the scaling factor for output of function values.
Definition: LineSearch.hpp:209
LineSearch lineSearch
We need a LineSearch object: backtracking line search.
#define IPOE(message)
Macro to allow file and line names in exceptions.
StoppingCondition
Enumeration for stopping conditions.
gsl::matrix modelHessian(gsl::matrix &matrix)
Find a Cholesky decomposition (both upper and lower triangles returned in matrix) of a matrix + mu * ...
Definition: GSL.cc:102
gsl::vector b
This vector is used in singular value decompostion and its size is set in constructor if it is needed...
StoppingCondition lastStoppingCondition
Last stopping condition.
ipo_function::DerivativesEstimates & de
The object containing the function and derivatives estimates.
Namespace for classes that handle details of interior-point optimisation that are not ordinarily acce...
Definition: Model.hpp:316
This class computes a function at a vector.
Definition: Function.hpp:38
void setGradTol(double const gradTol)
Set the value of gradTol.
gsl::vector S
This vector is used in singular value decompostion and its size is set in constructor if it is needed...
gsl::matrix V
This matrix is used in singular value decompostion and its size is set in constructor if it is needed...
gsl::matrix A
This matrix is used in singular value decompostion and its size is set in constructor if it is needed...
LineSearch::Parameters & getLineSearchParameters()
Get line search parameters by reference.
std::ostream * getOutputStream() const
Get the current value of outputStream.
Relative step size tolerance reached.
Line search parameters.
Definition: LineSearch.hpp:42
void setf_typ(double const f_typ)
Set the value of typical value of function.
double getGradTol() const
Get the value of gradTol.
virtual void setVector(gsl::vector const &vector)=0
Set the vector to a new value.
Parameters & getParameters()
Get parameters by reference.
This class allows us to test suitable stopping conditions.
void setPhaseI(bool const phaseI)
This function is designed to be used to help find an initial (phase I) feasible solution.