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