CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_REV_SPARSE_HES_INCLUDED 00003 # define CPPAD_REV_SPARSE_HES_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 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 RevSparseHes$$ 00018 $spell 00019 std 00020 VecAD 00021 Jacobian 00022 Jac 00023 Hessian 00024 Hes 00025 const 00026 Bool 00027 Dep 00028 proportional 00029 var 00030 $$ 00031 00032 $section Hessian Sparsity Pattern: Reverse Mode$$ 00033 00034 $index RevSparseHes$$ 00035 $index reverse, sparse Hessian$$ 00036 $index sparse, reverse Hessian$$ 00037 $index pattern, reverse Hessian$$ 00038 00039 $head Syntax$$ 00040 $icode%h% = %f%.RevSparseHes(%q%, %s%)%$$ 00041 00042 00043 $head Purpose$$ 00044 We use $latex F : \R^n \rightarrow \R^m$$ to denote the 00045 $xref/glossary/AD Function/AD function/$$ corresponding to $icode f$$. 00046 We define the matrix $latex H(x) \in \R^{q \times n}$$ 00047 as the partial with respect to $latex x$$, 00048 of the partial with respect to $latex u$$ (at $latex u = 0$$), 00049 of $latex S * F[ x + R * u ]$$ where 00050 $latex R \in \R^{n \times q}$$ and $latex S \in \R^{1 \times m}$$; i.e., 00051 $latex \[ 00052 H(x) = R^\T (S * F)^{(2)} ( x ) 00053 \] $$ 00054 Given a 00055 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00056 for the matrix $latex R$$ and the vector $latex S$$, 00057 $code RevSparseHes$$ returns a sparsity pattern for the $latex H(x)$$. 00058 00059 $head f$$ 00060 The object $icode f$$ has prototype 00061 $codei% 00062 const ADFun<%Base%> %f% 00063 %$$ 00064 00065 $head x$$ 00066 the sparsity pattern is valid for all values of the independent 00067 variables in $latex x \in \R^n$$ 00068 (even if it has $cref/CondExp/$$ or $cref/VecAD/$$ operations). 00069 00070 $head q$$ 00071 The argument $icode q$$ has prototype 00072 $codei% 00073 size_t %q% 00074 %$$ 00075 It specifies the number of columns in $latex R \in \R^{n \times q}$$ 00076 and the number of rows in 00077 $latex H(x) \in \R^{q \times n}$$. 00078 It must be the same value as in the previous $xref/ForSparseJac/$$ call 00079 $codei% 00080 %f%.ForSparseJac(%q%, %r%) 00081 %$$ 00082 00083 $head r$$ 00084 The argument $icode r$$ in the previous call 00085 $codei% 00086 %f%.ForSparseJac(%q%, %r%) 00087 %$$ 00088 is a $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00089 for the matrix $latex R$$ above. 00090 The type of the elements of 00091 $cref/VectorSet/RevSparseHes/VectorSet/$$ must be the 00092 same as the type of the elements of $icode r$$. 00093 00094 $head s$$ 00095 The argument $icode s$$ has prototype 00096 $codei% 00097 const %VectorSet%& %s% 00098 %$$ 00099 (see $xref/RevSparseHes/VectorSet/VectorSet/$$ below) 00100 If it has elements of type $code bool$$, 00101 its size is $latex m$$. 00102 If it has elements of type $code std::set<size_t>$$, 00103 its size is one and all the elements of $icode%s%[0]%$$ 00104 are between zero and $latex m - 1$$. 00105 It specifies a 00106 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00107 for the vector $icode S$$. 00108 00109 $head h$$ 00110 The result $icode h$$ has prototype 00111 $codei% 00112 %VectorSet%& %h% 00113 %$$ 00114 (see $xref/RevSparseHes/VectorSet/VectorSet/$$ below). 00115 If it has elements of type $code bool$$, 00116 its size is $latex q * n$$. 00117 If it has elements of type $code std::set<size_t>$$, 00118 its size is $latex q$$. 00119 It specifies a 00120 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ 00121 for the matrix $latex H(x)$$. 00122 00123 $head VectorSet$$ 00124 The type $icode VectorSet$$ must be a $cref/SimpleVector/$$ class with 00125 $xref/SimpleVector/Elements of Specified Type/elements of type/$$ 00126 $code bool$$ or $code std::set<size_t>$$; 00127 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00128 of the difference. 00129 The type of the elements of 00130 $cref/VectorSet/RevSparseHes/VectorSet/$$ must be the 00131 same as the type of the elements of $icode r$$. 00132 00133 $head Entire Sparsity Pattern$$ 00134 Suppose that $latex q = n$$ and 00135 $latex R \in \R^{n \times n}$$ is the $latex n \times n$$ identity matrix. 00136 Further suppose that the $latex S$$ is the $th k$$ 00137 $xref/glossary/Elementary Vector/elementary vector/$$; i.e. 00138 $latex \[ 00139 S_j = \left\{ \begin{array}{ll} 00140 1 & {\rm if} \; j = k 00141 \\ 00142 0 & {\rm otherwise} 00143 \end{array} \right. 00144 \] $$ 00145 In this case, 00146 the corresponding value $icode h$$ is a 00147 sparsity pattern for the Hessian matrix 00148 $latex F_k^{(2)} (x) \in \R^{n \times n}$$. 00149 00150 $head Example$$ 00151 $children% 00152 example/rev_sparse_hes.cpp 00153 %$$ 00154 The file 00155 $xref/RevSparseHes.cpp/$$ 00156 contains an example and test of this operation. 00157 It returns true if it succeeds and false otherwise. 00158 00159 $end 00160 ----------------------------------------------------------------------------- 00161 */ 00162 # include <algorithm> 00163 00164 CPPAD_BEGIN_NAMESPACE 00165 /*! 00166 \file rev_sparse_hes.hpp 00167 Reverse mode Hessian sparsity patterns. 00168 */ 00169 00170 /*! 00171 Calculate Hessian sparsity patterns using reverse mode. 00172 00173 The C++ source code corresponding to this operation is 00174 \verbatim 00175 h = f.RevSparseHes(q, s) 00176 \endverbatim 00177 00178 \tparam Base 00179 is the base type for this recording. 00180 00181 \tparam VectorSet 00182 is a simple vector with elements of type \c bool. 00183 00184 \tparam Sparsity 00185 is either \c sparse_pack or \c sparse_set. 00186 00187 \param q 00188 is the value of \a q in the 00189 by the previous call of the form 00190 \verbatim 00191 f.ForSparseJac(q, r) 00192 \endverbatim 00193 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00194 00195 \param s 00196 is a vector with size \c m that specifies the sparsity pattern 00197 for the vector \f$ S \f$, 00198 where \c m is the number of dependent variables 00199 corresponding to the operation sequence stored in \a play. 00200 00201 \param h 00202 the input value of \a h must be a vector with size \c q*n. 00203 The input value of its elements does not matter. 00204 On output, \a h is the sparsity pattern for the matrix 00205 \f[ 00206 H(x) = R^T (S * F)^{(2)} (x) 00207 \f] 00208 where \f$ F \f$ is the function corresponding to the operation sequence 00209 and \a x is any argument value. 00210 00211 \param total_num_var 00212 is the total number of variables in this recording. 00213 00214 \param dep_taddr 00215 maps dependendent variable index 00216 to the corresponding variable in the tape. 00217 00218 \param ind_taddr 00219 maps independent variable index 00220 to the corresponding variable in the tape. 00221 00222 \param play 00223 is the recording that defines the function we are computing the sparsity 00224 pattern for. 00225 00226 \param for_jac_sparsity 00227 is a vector of sets containing the 00228 the forward Jacobian sparsity pattern corresponding to 00229 $latex R$$ for all of the variables on the tape. 00230 */ 00231 00232 template <class Base, class VectorSet, class Sparsity> 00233 void RevSparseHesBool( 00234 size_t q , 00235 const VectorSet& s , 00236 VectorSet& h , 00237 size_t total_num_var , 00238 CppAD::vector<size_t>& dep_taddr , 00239 CppAD::vector<size_t>& ind_taddr , 00240 CppAD::player<Base>& play , 00241 Sparsity& for_jac_sparsity ) 00242 { 00243 // temporary indices 00244 size_t i, j; 00245 00246 // check Vector is Simple VectorSet class with bool elements 00247 CheckSimpleVector<bool, VectorSet>(); 00248 00249 // range and domain dimensions for F 00250 size_t m = dep_taddr.size(); 00251 size_t n = ind_taddr.size(); 00252 00253 CPPAD_ASSERT_KNOWN( 00254 q == for_jac_sparsity.end(), 00255 "RevSparseHes: q (first argument) is not equal to its value\n" 00256 "in the previous call to ForSparseJac with this ADFun object." 00257 ); 00258 CPPAD_ASSERT_KNOWN( 00259 s.size() == m, 00260 "RevSparseHes: s (second argument) length is not equal to\n" 00261 "range dimension for ADFun object." 00262 ); 00263 00264 // Array that will hold reverse Jacobian dependency flag. 00265 // Initialize as true for the dependent variables. 00266 bool *RevJac = CPPAD_NULL; 00267 RevJac = CPPAD_TRACK_NEW_VEC(total_num_var, RevJac); 00268 for(i = 0; i < total_num_var; i++) 00269 RevJac[i] = false; 00270 for(i = 0; i < m; i++) 00271 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00272 RevJac[ dep_taddr[i] ] = s[i]; 00273 } 00274 00275 00276 // vector of sets that will hold reverse Hessain values 00277 Sparsity rev_hes_sparsity; 00278 rev_hes_sparsity.resize(total_num_var, q); 00279 00280 // compute the Hessian sparsity patterns 00281 RevHesSweep( 00282 n, 00283 total_num_var, 00284 &play, 00285 for_jac_sparsity, 00286 RevJac, 00287 rev_hes_sparsity 00288 ); 00289 00290 // return values corresponding to independent variables 00291 CPPAD_ASSERT_UNKNOWN( h.size() == n * q ); 00292 00293 // j is index corresponding to reverse mode partial 00294 for(j = 0; j < n; j++) 00295 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var ); 00296 00297 // ind_taddr[j] is operator taddr for j-th independent variable 00298 CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == j + 1 ); 00299 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00300 00301 // extract the result from rev_hes_sparsity 00302 for(i = 0; i < q; i++) 00303 h[ i * n + j ] = false; 00304 CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q ); 00305 rev_hes_sparsity.begin(j + 1); 00306 i = rev_hes_sparsity.next_element(); 00307 while( i < q ) 00308 { h[ i * n + j ] = true; 00309 i = rev_hes_sparsity.next_element(); 00310 } 00311 } 00312 00313 // free local memory used for the calculation 00314 CPPAD_TRACK_DEL_VEC(RevJac); 00315 00316 return; 00317 } 00318 00319 /*! 00320 Calculate Hessian sparsity patterns using reverse mode. 00321 00322 The C++ source code corresponding to this operation is 00323 \verbatim 00324 h = f.RevSparseHes(q, s) 00325 \endverbatim 00326 00327 \tparam Base 00328 is the base type for this recording. 00329 00330 \tparam VectorSet 00331 is a simple vector with elements of type \c std::set<size_t>. 00332 00333 \tparam Sparsity 00334 is either \c sparse_pack or \c sparse_set. 00335 00336 \param q 00337 is the value of \a q in the 00338 by the previous call of the form 00339 \verbatim 00340 f.ForSparseJac(q, r) 00341 \endverbatim 00342 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00343 00344 \param s 00345 is a vector with size \c m that specifies the sparsity pattern 00346 for the vector \f$ S \f$, 00347 where \c m is the number of dependent variables 00348 corresponding to the operation sequence stored in \a play. 00349 00350 \param h 00351 the input value of \a h must be a vector with size \a q. 00352 On input, each element of \a h must be an empty set. 00353 The input value of its elements does not matter. 00354 On output, \a h is the sparsity pattern for the matrix 00355 \f[ 00356 H(x) = R^T (S * F)^{(2)} (x) 00357 \f] 00358 where \f$ F \f$ is the function corresponding to the operation sequence 00359 and \a x is any argument value. 00360 00361 \param total_num_var 00362 is the total number of variables in this recording. 00363 00364 \param dep_taddr 00365 maps dependendent variable index 00366 to the corresponding variable in the tape. 00367 00368 \param ind_taddr 00369 maps independent variable index 00370 to the corresponding variable in the tape. 00371 00372 \param play 00373 is the recording that defines the function we are computing the sparsity 00374 pattern for. 00375 00376 \param for_jac_sparsity 00377 is a vector of sets containing the 00378 the forward Jacobian sparsity pattern corresponding to 00379 $latex R$$ for all of the variables on the tape. 00380 */ 00381 00382 template <class Base, class VectorSet, class Sparsity> 00383 void RevSparseHesSet( 00384 size_t q , 00385 const VectorSet& s , 00386 VectorSet& h , 00387 size_t total_num_var , 00388 CppAD::vector<size_t>& dep_taddr , 00389 CppAD::vector<size_t>& ind_taddr , 00390 CppAD::player<Base>& play , 00391 Sparsity& for_jac_sparsity ) 00392 { 00393 // temporary indices 00394 size_t i, j; 00395 std::set<size_t>::const_iterator itr; 00396 00397 // check VectorSet is Simple Vector class with sets for elements 00398 static std::set<size_t> two, three; 00399 if( two.empty() ) 00400 { two.insert(2); 00401 three.insert(3); 00402 } 00403 CPPAD_ASSERT_UNKNOWN( two.size() == 1 ); 00404 CPPAD_ASSERT_UNKNOWN( three.size() == 1 ); 00405 CheckSimpleVector<std::set<size_t>, VectorSet>(two, three); 00406 00407 // range and domain dimensions for F 00408 # ifndef NDEBUG 00409 size_t m = dep_taddr.size(); 00410 # endif 00411 size_t n = ind_taddr.size(); 00412 00413 CPPAD_ASSERT_KNOWN( 00414 q == for_jac_sparsity.end(), 00415 "RevSparseHes: q (first argument) is not equal to its value\n" 00416 "in the previous call to ForSparseJac with this ADFun object." 00417 ); 00418 CPPAD_ASSERT_KNOWN( 00419 s.size() == 1, 00420 "RevSparseHes: s (second argument) length is not equal to one." 00421 ); 00422 00423 // Array that will hold reverse Jacobian dependency flag. 00424 // Initialize as true for the dependent variables. 00425 bool *RevJac = CPPAD_NULL; 00426 RevJac = CPPAD_TRACK_NEW_VEC(total_num_var, RevJac); 00427 for(i = 0; i < total_num_var; i++) 00428 RevJac[i] = false; 00429 itr = s[0].begin(); 00430 while( itr != s[0].end() ) 00431 { i = *itr++; 00432 CPPAD_ASSERT_KNOWN( 00433 i < m, 00434 "RevSparseHes: an element of the set s[0] has value " 00435 "greater than or equal m" 00436 ); 00437 CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00438 RevJac[ dep_taddr[i] ] = true; 00439 } 00440 00441 00442 // vector of sets that will hold reverse Hessain values 00443 Sparsity rev_hes_sparsity; 00444 rev_hes_sparsity.resize(total_num_var, q); 00445 00446 // compute the Hessian sparsity patterns 00447 RevHesSweep( 00448 n, 00449 total_num_var, 00450 &play, 00451 for_jac_sparsity, 00452 RevJac, 00453 rev_hes_sparsity 00454 ); 00455 00456 // return values corresponding to independent variables 00457 // j is index corresponding to reverse mode partial 00458 CPPAD_ASSERT_UNKNOWN( h.size() == q ); 00459 for(j = 0; j < n; j++) 00460 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var ); 00461 CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == j + 1 ); 00462 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00463 00464 // extract the result from rev_hes_sparsity 00465 // and add corresponding elements to result sets in h 00466 CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q ); 00467 rev_hes_sparsity.begin(j+1); 00468 i = rev_hes_sparsity.next_element(); 00469 while( i < q ) 00470 { h[i].insert(j); 00471 i = rev_hes_sparsity.next_element(); 00472 } 00473 } 00474 00475 // free local memory used for the calculation 00476 CPPAD_TRACK_DEL_VEC(RevJac); 00477 00478 return; 00479 } 00480 00481 00482 00483 /*! 00484 User API for Hessian sparsity patterns using reverse mode. 00485 00486 The C++ source code corresponding to this operation is 00487 \verbatim 00488 h = f.RevSparseHes(q, r) 00489 \endverbatim 00490 00491 \tparam Base 00492 is the base type for this recording. 00493 00494 \tparam VectorSet 00495 is a simple vector with elements of type \c bool 00496 or \c std::set<size_t>. 00497 00498 \param q 00499 is the value of \a q in the 00500 by the previous call of the form 00501 \verbatim 00502 f.ForSparseJac(q, r, packed) 00503 \endverbatim 00504 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00505 The type of the element of \c r for the previous call to \c ForSparseJac 00506 must be the same as the type of the elements of \c s. 00507 00508 \param s 00509 is a vector with size \c m that specifies the sparsity pattern 00510 for the vector \f$ S \f$, 00511 where \c m is the number of dependent variables 00512 corresponding to the operation sequence stored in \a play. 00513 00514 \return 00515 is a vector with size \c q*n. 00516 containing a sparsity pattern for the matrix 00517 \f[ 00518 H(x) = R^T ( S * F)^{(2)} (x) 00519 \f] 00520 where \f$ F \f$ is the function corresponding to the operation sequence 00521 and \a x is any argument value. 00522 */ 00523 00524 template <class Base> 00525 template <class VectorSet> 00526 VectorSet ADFun<Base>::RevSparseHes(size_t q, const VectorSet& s) 00527 { VectorSet h; 00528 typedef typename VectorSet::value_type Set_type; 00529 00530 RevSparseHesCase( 00531 Set_type() , 00532 q , 00533 s , 00534 h 00535 ); 00536 00537 return h; 00538 } 00539 00540 /*! 00541 Private helper function for RevSparseHes(q, s). 00542 00543 All of the description in the public member function RevSparseHes(q, s) 00544 applies. 00545 00546 \param set_type 00547 is a \c bool value. This argument is used to dispatch to the proper source 00548 code depending on the vlaue of \c VectorSet::value_type. 00549 00550 \param q 00551 See \c RevSparseHes(q, s). 00552 00553 \param s 00554 See \c RevSparseHes(q, s). 00555 00556 \param h 00557 is the return value for the corresponging call to \c RevSparseJac(q, s). 00558 */ 00559 template <class Base> 00560 template <class VectorSet> 00561 void ADFun<Base>::RevSparseHesCase( 00562 bool set_type , 00563 size_t q , 00564 const VectorSet& s , 00565 VectorSet& h ) 00566 { size_t n = Domain(); 00567 h.resize(q * n ); 00568 00569 CPPAD_ASSERT_KNOWN( 00570 for_jac_sparse_pack_.n_set() > 0, 00571 "RevSparseHes: previous stored call to ForSparseJac did not " 00572 "use bool for the elements of r." 00573 ); 00574 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 ); 00575 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == total_num_var_ ); 00576 00577 // use sparse_pack for the calculation 00578 CppAD::RevSparseHesBool( 00579 q , 00580 s , 00581 h , 00582 total_num_var_ , 00583 dep_taddr_ , 00584 ind_taddr_ , 00585 play_ , 00586 for_jac_sparse_pack_ 00587 ); 00588 } 00589 00590 /*! 00591 Private helper function for RevSparseHes(q, s). 00592 00593 All of the description in the public member function RevSparseHes(q, s) 00594 applies. 00595 00596 \param set_type 00597 is a \c std::set<size_t> value. 00598 This argument is used to dispatch to the proper source 00599 code depending on the vlaue of \c VectorSet::value_type. 00600 00601 \param q 00602 See \c RevSparseHes(q, s). 00603 00604 \param s 00605 See \c RevSparseHes(q, s). 00606 00607 \param h 00608 is the return value for the corresponging call to \c RevSparseJac(q, s). 00609 */ 00610 template <class Base> 00611 template <class VectorSet> 00612 void ADFun<Base>::RevSparseHesCase( 00613 const std::set<size_t>& set_type , 00614 size_t q , 00615 const VectorSet& s , 00616 VectorSet& h ) 00617 { h.resize(q); 00618 00619 CPPAD_ASSERT_KNOWN( 00620 for_jac_sparse_set_.n_set() > 0, 00621 "RevSparseHes: previous stored call to ForSparseJac did not " 00622 "use std::set<size_t> for the elements of r." 00623 ); 00624 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 ); 00625 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == total_num_var_ ); 00626 00627 // use sparse_pack for the calculation 00628 CppAD::RevSparseHesSet( 00629 q , 00630 s , 00631 h , 00632 total_num_var_ , 00633 dep_taddr_ , 00634 ind_taddr_ , 00635 play_ , 00636 for_jac_sparse_set_ 00637 ); 00638 } 00639 00640 CPPAD_END_NAMESPACE 00641 # endif