CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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