CppAD: A C++ Algorithmic Differentiation Package  20130102
hessian.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_HESSIAN_INCLUDED
00003 # define CPPAD_HESSIAN_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin Hessian$$
00018 $spell
00019      hes
00020      typename
00021      Taylor
00022      HesLuDet
00023      const
00024 $$
00025 
00026 $index Hessian, driver$$
00027 $index second, derivative$$
00028 $index driver, Hessian$$
00029 
00030 $section Hessian: Easy Driver$$
00031 
00032 $head Syntax$$
00033 $icode%hes% = %f%.Hessian(%x%, %w%)
00034 %$$
00035 $icode%hes% = %f%.Hessian(%x%, %l%)
00036 %$$
00037 
00038 
00039 $head Purpose$$
00040 We use $latex F : B^n \rightarrow B^m$$ to denote the
00041 $cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
00042 The syntax above sets $icode hes$$ to the Hessian 
00043 The syntax above sets $icode h$$ to the Hessian 
00044 $latex \[
00045      hes = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 
00046 \] $$
00047 The routine $cref sparse_hessian$$ may be faster in the case
00048 where the Hessian is sparse.
00049 
00050 $head f$$
00051 The object $icode f$$ has prototype
00052 $codei%
00053      ADFun<%Base%> %f%
00054 %$$
00055 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
00056 (see $cref/Hessian Uses Forward/Hessian/Hessian Uses Forward/$$ below).
00057 
00058 $head x$$
00059 The argument $icode x$$ has prototype
00060 $codei%
00061      const %Vector% &%x%
00062 %$$
00063 (see $cref/Vector/Hessian/Vector/$$ below)
00064 and its size 
00065 must be equal to $icode n$$, the dimension of the
00066 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00067 It specifies
00068 that point at which to evaluate the Hessian.
00069 
00070 $head l$$
00071 If the argument $icode l$$ is present, it has prototype
00072 $codei%
00073      size_t %l%
00074 %$$
00075 and is less than $icode m$$, the dimension of the
00076 $cref/range/seq_property/Range/$$ space for $icode f$$.
00077 It specifies the component of $icode F$$
00078 for which we are evaluating the Hessian.
00079 To be specific, in the case where the argument $icode l$$ is present,
00080 $latex \[
00081      w_i = \left\{ \begin{array}{ll}
00082           1 & i = l \\
00083           0 & {\rm otherwise}
00084      \end{array} \right.
00085 \] $$
00086 
00087 $head w$$
00088 If the argument $icode w$$ is present, it has prototype
00089 $codei%
00090      const %Vector% &%w%
00091 %$$
00092 and size $latex m$$.
00093 It specifies the value of $latex w_i$$ in the expression 
00094 for $icode h$$.
00095 
00096 $head hes$$
00097 The result $icode hes$$ has prototype
00098 $codei%
00099      %Vector% %hes%
00100 %$$
00101 (see $cref/Vector/Hessian/Vector/$$ below)
00102 and its size is $latex n * n$$.
00103 For $latex j = 0 , \ldots , n - 1 $$ 
00104 and $latex \ell = 0 , \ldots , n - 1$$
00105 $latex \[
00106      hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
00107 \] $$
00108 
00109 $head Vector$$
00110 The type $icode Vector$$ must be a $cref SimpleVector$$ class with
00111 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
00112 $icode Base$$.
00113 The routine $cref CheckSimpleVector$$ will generate an error message
00114 if this is not the case.
00115 
00116 $head Hessian Uses Forward$$
00117 After each call to $cref Forward$$,
00118 the object $icode f$$ contains the corresponding 
00119 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00120 After a call to $code Hessian$$,
00121 the zero order Taylor coefficients correspond to
00122 $icode%f%.Forward(0, %x%)%$$
00123 and the other coefficients are unspecified.
00124 
00125 $head Example$$
00126 $children%
00127      example/hessian.cpp%
00128      example/hes_lagrangian.cpp
00129 %$$
00130 The routines 
00131 $cref hessian.cpp$$ and
00132 $cref hes_lagrangian.cpp$$
00133 are examples and tests of $code Hessian$$.
00134 They return $code true$$, if they succeed and $code false$$ otherwise.
00135 
00136 
00137 $end
00138 -----------------------------------------------------------------------------
00139 */
00140 
00141 //  BEGIN CppAD namespace
00142 namespace CppAD {
00143 
00144 template <typename Base>
00145 template <typename Vector>
00146 Vector ADFun<Base>::Hessian(const Vector &x, size_t l)
00147 {    size_t i, m = Range();
00148      CPPAD_ASSERT_KNOWN(
00149           l < m,
00150           "Hessian: index i is not less than range dimension for f"
00151      );
00152 
00153      Vector w(m);
00154      for(i = 0; i < m; i++)
00155           w[i] = Base(0);
00156      w[l] = Base(1);
00157 
00158      return Hessian(x, w);
00159 }
00160 
00161 
00162 template <typename Base>
00163 template <typename Vector>
00164 Vector ADFun<Base>::Hessian(const Vector &x, const Vector &w)
00165 {    size_t j;
00166      size_t k;
00167 
00168      size_t n = Domain();
00169 
00170      // check Vector is Simple Vector class with Base type elements
00171      CheckSimpleVector<Base, Vector>();
00172 
00173      CPPAD_ASSERT_KNOWN(
00174           size_t(x.size()) == n,
00175           "Hessian: length of x not equal domain dimension for f"
00176      ); 
00177      CPPAD_ASSERT_KNOWN(
00178           size_t(w.size()) == Range(),
00179           "Hessian: length of w not equal range dimension for f"
00180      );
00181 
00182      // point at which we are evaluating the Hessian
00183      Forward(0, x);
00184 
00185      // define the return value
00186      Vector hes(n * n);
00187 
00188      // direction vector for calls to forward
00189      Vector u(n);
00190      for(j = 0; j < n; j++)
00191           u[j] = Base(0);
00192 
00193 
00194      // location for return values from Reverse
00195      Vector ddw(n * 2);
00196 
00197      // loop over forward directions
00198      for(j = 0; j < n; j++)
00199      {    // evaluate partials of entire function w.r.t. j-th coordinate
00200           u[j] = Base(1);
00201           Forward(1, u);
00202           u[j] = Base(0);
00203 
00204           // evaluate derivative of partial corresponding to F_i
00205           ddw = Reverse(2, w);
00206 
00207           // return desired components
00208           for(k = 0; k < n; k++)
00209                hes[k * n + j] = ddw[k * 2 + 1];
00210      }
00211           
00212      return hes;
00213 }
00214 
00215 } // END CppAD namespace
00216 
00217 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines