CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_HESSIAN_INCLUDED 00003 # define CPPAD_SPARSE_HESSIAN_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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 sparse_hessian$$ 00018 $spell 00019 CppAD 00020 valarray 00021 std 00022 Bool 00023 hes 00024 const 00025 Taylor 00026 $$ 00027 00028 $section Sparse Hessian: Easy Driver$$ 00029 $index SparseHessian$$ 00030 $index hessian, sparse$$ 00031 00032 $head Syntax$$ 00033 $codei%%hes% = %f%.SparseHessian(%x%, %w%) 00034 %$$ 00035 $codei%%hes% = %f%.SparseHessian(%x%, %w%, %p%)%$$ 00036 00037 $head Purpose$$ 00038 We use $latex F : \R^n \rightarrow \R^m$$ do denote the 00039 $cref/AD function/glossary/AD Function/$$ 00040 corresponding to $icode f$$. 00041 The syntax above sets $icode hes$$ to the Hessian 00042 $latex \[ 00043 hes = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 00044 \] $$ 00045 This routine assumes that the matrix $icode hes$$ is sparse 00046 and uses this assumption sparse to reduce the amount of computation necessary. 00047 One should use speed tests to verify that results are computed faster 00048 than when using the routine $cref/Hessian/$$. 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/Uses Forward/sparse_hessian/Uses Forward/$$ below). 00057 00058 $head x$$ 00059 The argument $icode x$$ has prototype 00060 $codei% 00061 const %VectorBase%& %x% 00062 %$$ 00063 (see $cref/VectorBase/sparse_hessian/VectorBase/$$ 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 w$$ 00071 The argument $icode w$$ has prototype 00072 $codei% 00073 const %VectorBase%& %w% 00074 %$$ 00075 and size $latex m$$. 00076 It specifies the value of $latex w_i$$ in the expression 00077 for $icode hes$$. 00078 The more components of $latex w$$ that are identically zero, 00079 the more sparse the resulting Hessian may be (and hence the more efficient 00080 the calculation of $icode hes$$ may be). 00081 00082 $head p$$ 00083 The argument $icode p$$ is optional and has prototype 00084 $codei% 00085 const %VectorSet%& %p% 00086 %$$ 00087 (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below) 00088 If it has elements of type $code bool$$, 00089 its size is $latex n * n$$. 00090 If it has elements of type $code std::set<size_t>$$, 00091 its size is $latex n$$ and all its set elements are between 00092 zero and $latex n - 1$$. 00093 It specifies a 00094 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00095 for the Hessian $icode hes$$. 00096 $pre 00097 00098 $$ 00099 If this sparsity pattern does not change between calls to 00100 $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and 00101 pass this argument to $codei SparseHessian$$. 00102 In addition, 00103 if you specify $icode p$$, CppAD will use the same 00104 type of sparsity representation 00105 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00106 for its internal calculations. 00107 Otherwise, the representation 00108 for the internal calculations is unspecified. 00109 00110 $head hes$$ 00111 The result $icode hes$$ has prototype 00112 $codei% 00113 %VectorBase% %hes% 00114 %$$ 00115 and its size is $latex n * n$$. 00116 For $latex j = 0 , \ldots , n - 1 $$ 00117 and $latex \ell = 0 , \ldots , n - 1$$ 00118 $latex \[ 00119 hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x ) 00120 \] $$ 00121 00122 $head VectorBase$$ 00123 The type $icode VectorBase$$ must be a $cref/SimpleVector/$$ class with 00124 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00125 $icode Base$$. 00126 The routine $cref/CheckSimpleVector/$$ will generate an error message 00127 if this is not the case. 00128 00129 $head VectorSet$$ 00130 The type $icode VectorSet$$ must be a $xref/SimpleVector/$$ class with 00131 $xref/SimpleVector/Elements of Specified Type/elements of type/$$ 00132 $code bool$$ or $code std::set<size_t>$$; 00133 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00134 of the difference. 00135 The routine $xref/CheckSimpleVector/$$ will generate an error message 00136 if this is not the case. 00137 00138 $subhead Restrictions$$ 00139 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00140 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00141 corresponding set. 00142 According to section 26.3.2.3 of the 1998 C++ standard, 00143 $code std::valarray< std::set<size_t> >$$ does not satisfy 00144 this condition. 00145 00146 $head Uses Forward$$ 00147 After each call to $cref/Forward/$$, 00148 the object $icode f$$ contains the corresponding 00149 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00150 After $code SparseHessian$$, 00151 the previous calls to $xref/Forward/$$ are undefined. 00152 00153 $head Example$$ 00154 $children% 00155 example/sparse_hessian.cpp 00156 %$$ 00157 The routine 00158 $cref/sparse_hessian.cpp/$$ 00159 is examples and tests of $code sparse_hessian$$. 00160 It return $code true$$, if it succeeds and $code false$$ otherwise. 00161 00162 00163 $end 00164 ----------------------------------------------------------------------------- 00165 */ 00166 00167 CPPAD_BEGIN_NAMESPACE 00168 /*! 00169 \file sparse_hessian.hpp 00170 Sparse Hessian driver routine and helper functions. 00171 */ 00172 00173 /*! 00174 Private helper function for SparseHessian(x, w, p). 00175 00176 All of the description in the public member function SparseHessian(x, w, p) 00177 applies. 00178 00179 \param set_type 00180 is a \c bool value. This argument is used to dispatch to the proper source 00181 code depending on the value of \c VectorSet::value_type. 00182 00183 \param x 00184 See \c SparseHessian(x, w, p). 00185 00186 \param w 00187 See \c SparseHessian(x, w, p). 00188 00189 \param p 00190 See \c SparseHessian(x, w, p). 00191 00192 \param hes 00193 is the return value for the corresponding call to \c SparseHessian(x, w, p). 00194 On input, it must have size equal to the domain times range dimension 00195 for this ADFun<Base> object. 00196 On return, it will contain the Hessian. 00197 */ 00198 00199 template <class Base> 00200 template <class VectorBase, class VectorSet> 00201 void ADFun<Base>::SparseHessianCase( 00202 bool set_type , 00203 const VectorBase& x , 00204 const VectorBase& w , 00205 const VectorSet& p , 00206 VectorBase& hes ) 00207 { 00208 typedef CppAD::vector<size_t> SizeVector; 00209 typedef CppAD::vectorBool VectorBool; 00210 size_t j, k, l; 00211 00212 size_t n = Domain(); 00213 00214 // check Vector is Simple Vector class with bool elements 00215 CheckSimpleVector<bool, VectorSet>(); 00216 00217 // check Vector is Simple Vector class with Base type elements 00218 CheckSimpleVector<Base, VectorBase>(); 00219 00220 CPPAD_ASSERT_KNOWN( 00221 x.size() == n, 00222 "SparseHessian: size of x not equal domain dimension for f" 00223 ); 00224 CPPAD_ASSERT_KNOWN( 00225 w.size() == Range(), 00226 "SparseHessian: size of w not equal range dimension for f" 00227 ); 00228 CPPAD_ASSERT_KNOWN( 00229 p.size() == n * n, 00230 "SparseHessian: using bool values and size of p " 00231 "not equal square of domain dimension for f" 00232 ); 00233 00234 // initial coloring 00235 SizeVector color(n); 00236 for(j = 0; j < n; j++) 00237 color[j] = j; 00238 00239 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00240 // Graph Coloring in Optimization Revisited by 00241 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00242 VectorBool forbidden(n); 00243 for(j = 0; j < n; j++) 00244 { // initial all colors as ok for this column 00245 for(k = 0; k < n; k++) 00246 forbidden[k] = false; 00247 // for each row that is connected to column j 00248 for(l = 0; l < n; l++) if( p[l * n + j] ) 00249 { // for each column that is connected to row l 00250 for(k = 0; k < n; k++) if( p[l * n + k] & (j != k) ) 00251 forbidden[ color[k] ] = true; 00252 } 00253 k = 0; 00254 while( forbidden[k] && k < n ) 00255 { k++; 00256 CPPAD_ASSERT_UNKNOWN( k < n ); 00257 } 00258 color[j] = k; 00259 } 00260 size_t n_color = 1; 00261 for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1); 00262 00263 // some values 00264 const Base zero(0); 00265 const Base one(1); 00266 00267 // point at which we are evaluating the Hessian 00268 Forward(0, x); 00269 00270 // initialize the return value 00271 for(j = 0; j < n; j++) 00272 for(k = 0; k < n; k++) 00273 hes[j * n + k] = zero; 00274 00275 // direction vector for calls to forward 00276 VectorBase u(n); 00277 00278 // location for return values from Reverse 00279 VectorBase ddw(n * 2); 00280 00281 // loop over colors 00282 size_t c; 00283 for(c = 0; c < n_color; c++) 00284 { // determine all the colums with this color 00285 for(j = 0; j < n; j++) 00286 { if( color[j] == c ) 00287 u[j] = one; 00288 else u[j] = zero; 00289 } 00290 // call forward mode for all these columns at once 00291 Forward(1, u); 00292 00293 // evaluate derivative of w^T * F'(x) * u 00294 ddw = Reverse(2, w); 00295 00296 // set the corresponding components of the result 00297 for(j = 0; j < n; j++) if( color[j] == c ) 00298 { for(l = 0; l < n; l++) 00299 if( p[ l * n + j ] ) 00300 hes[l * n + j] = ddw[l * 2 + 1]; 00301 } 00302 } 00303 } 00304 /*! 00305 Private helper function for SparseHessian(x, w, p). 00306 00307 All of the description in the public member function SparseHessian(x, w, p) 00308 applies. 00309 00310 \param set_type 00311 is a \c std::set<size_t> value. 00312 This argument is used to dispatch to the proper source 00313 code depending on the value of \c VectorSet::value_type. 00314 00315 \param x 00316 See \c SparseHessian(x, w, p). 00317 00318 \param w 00319 See \c SparseHessian(x, w, p). 00320 00321 \param p 00322 See \c SparseHessian(x, w, p). 00323 00324 \param hes 00325 is the return value for the corresponding call to \c SparseHessian(x, w, p). 00326 On input, it must have size equal to the domain times range dimension 00327 for this ADFun<Base> object. 00328 On return, it will contain the Hessian. 00329 */ 00330 00331 template <class Base> 00332 template <class VectorBase, class VectorSet> 00333 void ADFun<Base>::SparseHessianCase( 00334 const std::set<size_t>& set_type , 00335 const VectorBase& x , 00336 const VectorBase& w , 00337 const VectorSet& p , 00338 VectorBase& hes ) 00339 { 00340 typedef CppAD::vector<size_t> SizeVector; 00341 typedef CppAD::vectorBool VectorBool; 00342 size_t j, k, l; 00343 00344 size_t n = Domain(); 00345 00346 // check VectorSet is Simple Vector class with sets for elements 00347 static std::set<size_t> two, three; 00348 if( two.empty() ) 00349 { two.insert(2); 00350 three.insert(3); 00351 } 00352 CPPAD_ASSERT_UNKNOWN( two.size() == 1 ); 00353 CPPAD_ASSERT_UNKNOWN( three.size() == 1 ); 00354 CheckSimpleVector<std::set<size_t>, VectorSet>(two, three); 00355 00356 // check Vector is Simple Vector class with Base type elements 00357 CheckSimpleVector<Base, VectorBase>(); 00358 00359 CPPAD_ASSERT_KNOWN( 00360 x.size() == n, 00361 "SparseHessian: size of x not equal domain dimension for f" 00362 ); 00363 CPPAD_ASSERT_KNOWN( 00364 w.size() == Range(), 00365 "SparseHessian: size of w not equal range dimension for f" 00366 ); 00367 CPPAD_ASSERT_KNOWN( 00368 p.size() == n, 00369 "SparseHessian: using set values and size of p " 00370 "not equal domain dimension for f" 00371 ); 00372 00373 // initial coloring 00374 SizeVector color(n); 00375 for(j = 0; j < n; j++) 00376 color[j] = j; 00377 00378 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00379 // Graph Coloring in Optimization Revisited by 00380 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00381 VectorBool forbidden(n); 00382 std::set<size_t>::const_iterator itr_k, itr_l; 00383 for(j = 0; j < n; j++) 00384 { // initial all colors as ok for this column 00385 for(k = 0; k < n; k++) 00386 forbidden[k] = false; 00387 00388 // for each row that is connected to column j 00389 itr_l = p[j].begin(); 00390 while( itr_l != p[j].end() ) 00391 { l = *itr_l++; 00392 // for each column that is connected to row l 00393 // (same as each row that is connect to column l) 00394 itr_k = p[l].begin(); 00395 while( itr_k != p[l].end() ) 00396 { k = *itr_k++; 00397 forbidden[ color[k] ] = (j != k); 00398 } 00399 } 00400 k = 0; 00401 while( forbidden[k] && k < n ) 00402 { k++; 00403 CPPAD_ASSERT_UNKNOWN( k < n ); 00404 } 00405 color[j] = k; 00406 } 00407 size_t n_color = 1; 00408 for(k = 0; k < n; k++) n_color = std::max(n_color, color[k] + 1); 00409 00410 // some values 00411 const Base zero(0); 00412 const Base one(1); 00413 00414 // point at which we are evaluating the Hessian 00415 Forward(0, x); 00416 00417 // initialize the return value 00418 for(j = 0; j < n; j++) 00419 for(k = 0; k < n; k++) 00420 hes[j * n + k] = zero; 00421 00422 // direction vector for calls to forward 00423 VectorBase u(n); 00424 00425 // location for return values from Reverse 00426 VectorBase ddw(n * 2); 00427 00428 // loop over colors 00429 size_t c; 00430 for(c = 0; c < n_color; c++) 00431 { // determine all the colums with this color 00432 for(j = 0; j < n; j++) 00433 { if( color[j] == c ) 00434 u[j] = one; 00435 else u[j] = zero; 00436 } 00437 // call forward mode for all these columns at once 00438 Forward(1, u); 00439 00440 // evaluate derivative of w^T * F'(x) * u 00441 ddw = Reverse(2, w); 00442 00443 // set the corresponding components of the result 00444 for(j = 0; j < n; j++) if( color[j] == c ) 00445 { itr_l = p[j].begin(); 00446 while( itr_l != p[j].end() ) 00447 { l = *itr_l++; 00448 hes[l * n + j] = ddw[l * 2 + 1]; 00449 } 00450 } 00451 } 00452 } 00453 00454 /*! 00455 Compute a sparse Hessian. 00456 00457 The C++ source code coresponding to this operation is 00458 \verbatim 00459 hes = SparseHessian(x, w) 00460 \endverbatim 00461 00462 00463 \tparam Base 00464 is the base type for the recording that is stored in this 00465 ADFun<Base object. 00466 00467 \tparam VectorBase 00468 is a simple vector class with elements of the \a Base. 00469 00470 \tparam VectorSet 00471 is a simple vector class with elements of type 00472 \c bool or \c std::set<size_t>. 00473 00474 \param x 00475 is a vector specifing the point at which to compute the Hessian. 00476 00477 \param w 00478 The Hessian is computed for a weighted sum of the components 00479 of the function corresponding to this ADFun<Base> object. 00480 The argument \a w specifies the weights for each component. 00481 It must have size equal to the range dimension for this ADFun<Base> object. 00482 00483 \param p 00484 is a sparsity pattern for the Hessian. 00485 00486 \return 00487 Will be a vector of size \c n * n containing the Hessian of 00488 at the point specified by \a x 00489 (where \c n is the domain dimension for this ADFun<Base> object). 00490 */ 00491 template <class Base> 00492 template <class VectorBase, class VectorSet> 00493 VectorBase ADFun<Base>::SparseHessian( 00494 const VectorBase& x, const VectorBase& w, const VectorSet& p 00495 ) 00496 { size_t n = Domain(); 00497 VectorBase hes(n * n); 00498 00499 typedef typename VectorSet::value_type Set_type; 00500 SparseHessianCase(Set_type(), x, w, p, hes); 00501 00502 return hes; 00503 } 00504 00505 /*! 00506 Compute a sparse Hessian 00507 00508 The C++ source code coresponding to this operation is 00509 \verbatim 00510 hes = SparseHessian(x, w) 00511 \endverbatim 00512 00513 00514 \tparam Base 00515 is the base type for the recording that is stored in this 00516 ADFun<Base object. 00517 00518 \tparam VectorBase 00519 is a simple vector class with elements of the \a Base. 00520 00521 \param x 00522 is a vector specifing the point at which to compute the Hessian. 00523 00524 \param w 00525 The Hessian is computed for a weighted sum of the components 00526 of the function corresponding to this ADFun<Base> object. 00527 The argument \a w specifies the weights for each component. 00528 It must have size equal to the range dimension for this ADFun<Base> object. 00529 00530 \return 00531 Will be a vector of size \c n * n containing the Hessian of 00532 at the point specified by \a x 00533 (where \c n is the domain dimension for this ADFun<Base> object). 00534 */ 00535 template <class Base> 00536 template <class VectorBase> 00537 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w) 00538 { size_t i, j, k; 00539 typedef CppAD::vectorBool VectorBool; 00540 00541 size_t m = Range(); 00542 size_t n = Domain(); 00543 00544 // determine the sparsity pattern p for Hessian of w^T F 00545 VectorBool r(n * n); 00546 for(j = 0; j < n; j++) 00547 { for(k = 0; k < n; k++) 00548 r[j * n + k] = false; 00549 r[j * n + j] = true; 00550 } 00551 ForSparseJac(n, r); 00552 // 00553 VectorBool s(m); 00554 for(i = 0; i < m; i++) 00555 s[i] = w[i] != 0; 00556 VectorBool p = RevSparseHes(n, s); 00557 00558 // compute sparse Hessian 00559 VectorBase hes(n * n); 00560 bool set_type = true; // only used to dispatch compiler to proper case 00561 SparseHessianCase(set_type, x, w, p, hes); 00562 00563 return hes; 00564 } 00565 00566 CPPAD_END_NAMESPACE 00567 # endif