CppAD: A C++ Algorithmic Differentiation Package 20110419
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-09 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common 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 $syntax%%hes% = %f%.Hessian(%x%, %w%)
00034 %$$
00035 $syntax%%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 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$.
00042 The syntax above sets $italic hes$$ to the Hessian 
00043 The syntax above sets $italic 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 $italic f$$ has prototype
00052 $syntax%
00053         ADFun<%Base%> %f%
00054 %$$
00055 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$
00056 (see $xref/Hessian/Hessian Uses Forward/Hessian Uses Forward/$$ below).
00057 
00058 $head x$$
00059 The argument $italic x$$ has prototype
00060 $syntax%
00061         const %Vector% &%x%
00062 %$$
00063 (see $xref/Hessian/Vector/Vector/$$ below)
00064 and its size 
00065 must be equal to $italic n$$, the dimension of the
00066 $xref/seq_property/Domain/domain/$$ space for $italic f$$.
00067 It specifies
00068 that point at which to evaluate the Hessian.
00069 
00070 $head l$$
00071 If the argument $italic l$$ is present, it has prototype
00072 $syntax%
00073         size_t %l%
00074 %$$
00075 and is less than $italic m$$, the dimension of the
00076 $xref/seq_property/Range/range/$$ space for $italic f$$.
00077 It specifies the component of $italic F$$
00078 for which we are evaluating the Hessian.
00079 To be specific, in the case where the argument $italic 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 $italic w$$ is present, it has prototype
00089 $syntax%
00090         const %Vector% &%w%
00091 %$$
00092 and size $latex m$$.
00093 It specifies the value of $latex w_i$$ in the expression 
00094 for $italic h$$.
00095 
00096 $head hes$$
00097 The result $italic hes$$ has prototype
00098 $syntax%
00099         %Vector% %hes%
00100 %$$
00101 (see $xref/Hessian/Vector/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 $italic Vector$$ must be a $xref/SimpleVector/$$ class with
00111 $xref/SimpleVector/Elements of Specified Type/elements of type/$$
00112 $italic Base$$.
00113 The routine $xref/CheckSimpleVector/$$ will generate an error message
00114 if this is not the case.
00115 
00116 $head Hessian Uses Forward$$
00117 After each call to $xref/Forward/$$,
00118 the object $italic f$$ contains the corresponding 
00119 $xref/glossary/Taylor Coefficient/Taylor coefficients/$$.
00120 After $code Hessian$$,
00121 the previous calls to $xref/Forward/$$ are undefined.
00122 
00123 $head Example$$
00124 $children%
00125         example/hessian.cpp%
00126         example/hes_lagrangian.cpp
00127 %$$
00128 The routines 
00129 $cref/Hessian.cpp/$$ and
00130 $cref/HesLagrangian.cpp/$$
00131 are examples and tests of $code Hessian$$.
00132 They return $code true$$, if they succeed and $code false$$ otherwise.
00133 
00134 
00135 $end
00136 -----------------------------------------------------------------------------
00137 */
00138 
00139 //  BEGIN CppAD namespace
00140 namespace CppAD {
00141 
00142 template <typename Base>
00143 template <typename Vector>
00144 Vector ADFun<Base>::Hessian(const Vector &x, size_t l)
00145 {       size_t i, m = Range();
00146         CPPAD_ASSERT_KNOWN(
00147                 l < m,
00148                 "Hessian: index i is not less than range dimension for f"
00149         );
00150 
00151         Vector w(m);
00152         for(i = 0; i < m; i++)
00153                 w[i] = Base(0);
00154         w[l] = Base(1);
00155 
00156         return Hessian(x, w);
00157 }
00158 
00159 
00160 template <typename Base>
00161 template <typename Vector>
00162 Vector ADFun<Base>::Hessian(const Vector &x, const Vector &w)
00163 {       size_t j;
00164         size_t k;
00165 
00166         size_t n = Domain();
00167 
00168         // check Vector is Simple Vector class with Base type elements
00169         CheckSimpleVector<Base, Vector>();
00170 
00171         CPPAD_ASSERT_KNOWN(
00172                 x.size() == n,
00173                 "Hessian: length of x not equal domain dimension for f"
00174         ); 
00175         CPPAD_ASSERT_KNOWN(
00176                 w.size() == Range(),
00177                 "Hessian: length of w not equal range dimension for f"
00178         );
00179 
00180         // point at which we are evaluating the Hessian
00181         Forward(0, x);
00182 
00183         // define the return value
00184         Vector hes(n * n);
00185 
00186         // direction vector for calls to forward
00187         Vector u(n);
00188         for(j = 0; j < n; j++)
00189                 u[j] = Base(0);
00190 
00191 
00192         // location for return values from Reverse
00193         Vector ddw(n * 2);
00194 
00195         // loop over forward directions
00196         for(j = 0; j < n; j++)
00197         {       // evaluate partials of entire function w.r.t. j-th coordinate
00198                 u[j] = Base(1);
00199                 Forward(1, u);
00200                 u[j] = Base(0);
00201 
00202                 // evaluate derivative of partial corresponding to F_i
00203                 ddw = Reverse(2, w);
00204 
00205                 // return desired components
00206                 for(k = 0; k < n; k++)
00207                         hes[k * n + j] = ddw[k * 2 + 1];
00208         }
00209                 
00210         return hes;
00211 }
00212 
00213 } // END CppAD namespace
00214 
00215 # endif