29 double constexpr gradTolDefault
30 { std::exp( std::log( std::numeric_limits<double>::epsilon() ) / 3 ) };
31 double constexpr stepTolDefault { gradTolDefault * gradTolDefault };
54 gsl::vector
const equalityConstraintsBounds )
55 : equalityConstraintsMatrix { equalityConstraintsMatrix },
56 function( function ), de( de ), lineSearch {
function }
62 if( equalityConstraintsMatrix ){
63 size_t const SIZE1 { equalityConstraintsMatrix.size1() };
64 size_t const SIZE2 { equalityConstraintsMatrix.size2() };
65 if( SIZE1 > 0 and SIZE2 > 0 ){
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." );
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];
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." );
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." );
97 size_t const SIZE { x.size() };
98 gsl::vector delta( SIZE );
106 size_t iteration { 0 };
107 objectiveValue =
function( x );
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++ <<
" ";
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 ) <<
" ";
119 *outputStream <<
")";
121 size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
122 *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
127 double lastEntry { *x.rbegin() };
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 ) };
146 if( j > i )
A.set( j, i, d );
148 for(
size_t j { 0 }; j < SIZE1; ++j ){
150 A.set( SIZE2 + j, i, d );
151 A.set( i, SIZE2 + j, d );
154 for(
size_t i { SIZE2 }; i < SIZE1 + SIZE2; ++i )
155 for(
size_t j { i }; j < SIZE1 + SIZE2; ++j ){
157 if( j > i )
A.set( j, i, 0 );
160 auto handler = gsl::exception::set_handler_gsl_exceptions();
162 gsl::linalg::SV_decomp(
A,
V,
S,
work );
165 for(
size_t i { 0 }; i < SIZE2; ++i )
b[i] = g[i];
166 for(
size_t i { SIZE2 }; i < SIZE1 + SIZE2; ++i )
b[i] = 0;
167 gsl::linalg::SV_solve(
A,
V,
S,
b,
work );
169 for(
size_t i { 0 }; i < SIZE2; ++i ) delta[i] =
work[i];
170 }
catch( gsl::exception& e ){
174 gsl::exception::set_handler( handler );
175 throw IPOE(
"ipo::detail::NewtonDescent::operator(): "
176 "unknown exception encountered." );
178 gsl::exception::set_handler( handler );
182 gsl::linalg::cholesky_solve( choleDecomp,
de.
gradient(), delta );
189 objectiveValue =
function( x );
192 if( outputStream !=
nullptr ){
193 *outputStream <<
" " << std::setw( 4 ) << iteration <<
" ";
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 ) <<
" ";
200 *outputStream <<
")";
202 size_t const d { std::max( static_cast<size_t>( 6 ), digits ) };
203 *outputStream << std::setw( 4 + d ) << std::setprecision( d ) << std::fixed
207 double const xLastEntry { *x.rbegin() };
210 if( ls < std::numeric_limits<double>::epsilon() )
break;
213 if( xLastEntry < lastEntry ) lastEntry = xLastEntry;
220 if( outputStream !=
nullptr )
221 *outputStream <<
"Newton descent finished." << std::endl << std::endl;
222 return objectiveValue;
227 this->outputStream = outputStream;
236 : newtonDescent( newtonDescent ){}
241 gsl::vector
const& gradient, gsl::vector
const& vector,
242 gsl::vector
const& lastVector )
const {
244 gsl::vector
const& x_typ
252 if( iteration > const_cast<NewtonDescent&>( newtonDescent ).
getParameters().getIterMax() ){
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;
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 ) ) };
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() )
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.
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.
Gradient tolerance reached.
void setIterMax(size_t const iterMax)
Set the value of iterMax.
Parameters & getParameters()
Get parameters by reference.
void setFunctionScale(double const functionScale)
Set the scaling factor for output of function values.
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 * ...
gsl::vector b
This vector is used in singular value decompostion and its size is set in constructor if it is needed...
Maximum iterations reached.
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...
This class computes a function at a vector.
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.
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.