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