Optimo  0.1.0
A C++ header library for optimization
 All Classes Functions Variables Pages
numerical_hessian.h
1 // Copyright (C) 2013 Victor Fragoso <vfragoso@cs.ucsb.edu>
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted provided that the following conditions are
6 // met:
7 //
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 //
11 // * Redistributions in binary form must reproduce the above
12 // copyright notice, this list of conditions and the following
13 // disclaimer in the documentation and/or other materials provided
14 // with the distribution.
15 //
16 // * Neither the name of the University of California, Santa Barbara nor the
17 // names of its contributors may be used to endorse or promote products
18 // derived from this software without specific prior written permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23 // ARE DISCLAIMED. IN NO EVENT SHALL VICTOR FRAGOSO BE LIABLE FOR ANY DIRECT,
24 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
25 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
26 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
27 // ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 //
31 
32 #ifndef OPTIMO_CORE_NUMERICAL_HESSIAN_H_
33 #define OPTIMO_CORE_NUMERICAL_HESSIAN_H_
34 
35 #include <Eigen/Core>
36 #include "optimo/core/objects.h"
37 
38 namespace optimo {
40 
42 template <typename Scalar, uint n>
43 class NumericalHessian : public HessianFunctor<Scalar, n> {
44  public:
47  const Scalar h) :
48  HessianFunctor<Scalar, n>(), objective_(obj), h_(h) { }
49 
50  virtual ~NumericalHessian(void) { }
51 
53  virtual void
54  operator()(const Eigen::Matrix<Scalar, n, 1>& x,
55  Eigen::Matrix<Scalar, n, n>* H) const {
56  Scalar den = h_*h_;
57  // 1. Calculate Diagonal of the Hessian
58  Eigen::Matrix<Scalar, n, 1> xnn, xpp, xnp, xpn;
59  for (unsigned int i = 0; i < n; i++) {
60  xnn = x;
61  xpp = x;
62  xpp(i, 0) += h_;
63  xnn(i, 0) -= h_;
64  (*H)(i, i) = (objective_(xpp) - 2*objective_(x) + objective_(xnn))/den;
65  }
66 
67  // 2. Compute the upper triangular part of the Hessian
68  den *= 4;
69  for (unsigned int i = 0; i < n; i++) {
70  for (unsigned int j = i + 1; j < n; j++) {
71  xpp = x;
72  xnp = x;
73  xpn = x;
74  xnn = x;
75 
76  xpp(i, 0) += h_;
77  xpp(j, 0) += h_;
78 
79  xpn(i, 0) += h_;
80  xpn(j, 0) -= h_;
81 
82  xnp(i, 0) -= h_;
83  xnp(j, 0) += h_;
84 
85  xnn(i, 0) -= h_;
86  xnn(j, 0) -= h_;
87 
88  (*H)(i, j) = (objective_(xpp) - objective_(xpn) -
89  objective_(xnp) + objective_(xnn))/den;
90  (*H)(j, i) = (*H)(i, j);
91  }
92  }
93  }
94 
95  protected:
96  const ObjectiveFunctor<Scalar, n>& objective_;
97  const Scalar h_;
98 };
99 } // optimo
100 #endif // OPTIMO_CORE_NUMERICAL_HESSIAN_H_
Numerical Hessian computation.
Definition: numerical_hessian.h:43
Objective functor.
Definition: objects.h:54
NumericalHessian(const ObjectiveFunctor< Scalar, n > &obj, const Scalar h)
Constructor.
Definition: numerical_hessian.h:46
Dense Hessian functor.
Definition: objects.h:104
virtual void operator()(const Eigen::Matrix< Scalar, n, 1 > &x, Eigen::Matrix< Scalar, n, n > *H) const
Computes Hessian matrix.
Definition: numerical_hessian.h:54