Interior-point-optimisation  1.0-1
Interior-pointoptimisationlibrary
ForwardDifferenceDerivativesEstimates.cc
Go to the documentation of this file.
1 /*
2  * $Id: ForwardDifferenceDerivativesEstimates.cc 137 2013-06-29 15:10:31Z jdl3 $
3  * Copyright (C) 2011, 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 
25 
26 using namespace ipo_function::detail;
27 
29 ForwardDifferenceDerivativesEstimates( Function& function, double const h, double const h2 )
30  : ForwardDifferenceGradientEstimate { function, h }, h2 { h2 }{}
31 
32 void
34  // Call parent class
36  // Now find SIZE
37  const size_t SIZE { vector.size() };
38  // // Set up storage for f(x + h2 e_i) and f(x + h e_i + h e_j)
39  auto diffVector2 // For Hessian
40  = gsl::vector( SIZE );
41  auto diffMatrix = gsl::matrix( SIZE, SIZE );
42  // Fill diffVector2
43  for( size_t i { 0 }; i < SIZE; ++i ){
44  double tmp { vector.get( i ) };
45  const_cast<gsl::vector&>( vector ).set( i, tmp + h2 );
46  diffVector2.set( i, function( vector ) );
47  const_cast<gsl::vector&>( vector ).set( i, tmp );
48  }
49  // Fill diffMatrix: upper triangle only
50  for( size_t i { 0 }; i < SIZE; ++i ){
51  double tmp1 { vector.get( i ) };
52  const_cast<gsl::vector&>( vector ).set( i, tmp1 + 2 * h2 );
53  diffMatrix.set( i, i, function( vector ) );
54  const_cast<gsl::vector&>( vector ).set( i, tmp1 + h2 );
55  for( size_t j { i + 1 }; j < SIZE; ++j ){
56  double tmp2 { vector.get( j ) };
57  const_cast<gsl::vector&>( vector ).set( j, tmp2 + h2 );
58  diffMatrix.set( i, j, function( vector ) );
59  const_cast<gsl::vector&>( vector ).set( j, tmp2 );
60  }
61  const_cast<gsl::vector&>( vector ).set( i, tmp1 );
62  }
63 
64  // set up functionHessian
65  if( functionHessian.get() == nullptr or functionHessian.size1() != SIZE or
66  functionHessian.size2() != SIZE )
67  functionHessian = gsl::matrix( SIZE, SIZE );
68 
69  // calculate Hessian: forward difference quotient
70  double const H { h2 * h2 };
71  for( size_t i { 0 }; i < SIZE; ++i ){
72  double const diff { (diffMatrix.get( i, i ) - 2 * diffVector2.get( i )
73  + functionValue) / H };
74  functionHessian.set( i, i, diff );
75  for( size_t j { i + 1 }; j < SIZE; ++j ){
76  double const diff { (diffMatrix.get( i, j ) - diffVector2.get( i )
77  - diffVector2.get( j ) + functionValue ) / H };
78  functionHessian.set( i, j, diff );
79  functionHessian.set( j, i, diff );
80  }
81  }
82 }
virtual void setVector(gsl::vector const &vector)
Set the vector to a new value.
double functionValue
The function value.
virtual void setVector(gsl::vector const &vector)
Set the vector to a new value.
This class estimates a function value, gradient and Hessian at a given vector.
Namespace for details of ipo_function that are not normally needed to construct and solve a convex op...
gsl::matrix functionHessian
The Hessian value.
This class computes a function at a vector.
Definition: Function.hpp:38
double const h2
The quotient size for hessian estimates: default is fourth root of machine epsilon.
ForwardDifferenceDerivativesEstimates(Function &function, double const h=std::sqrt(std::numeric_limits< double >::epsilon()), double const h2=std::sqrt(std::sqrt(std::numeric_limits< double >::epsilon())))
Find forward difference quotient estimates.