CppAD: A C++ Algorithmic Differentiation Package
20130102
|
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-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 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 : \B{R}^n \rightarrow \B{R}^m$$ to denote the 00045 $cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$. 00046 We define the matrix $latex H(x) \in \B{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 \B{R}^{n \times q}$$ and $latex S \in \B{R}^{1 \times m}$$; i.e., 00051 $latex \[ 00052 H(x) = R^\R{T} (S * F)^{(2)} ( x ) 00053 \] $$ 00054 Given a 00055 $cref/sparsity pattern/glossary/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 \B{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 \B{R}^{n \times q}$$ 00076 and the number of rows in 00077 $latex H(x) \in \B{R}^{q \times n}$$. 00078 It must be the same value as in the previous $cref 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 $cref/sparsity pattern/glossary/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 $cref/VectorSet/RevSparseHes/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 $cref/sparsity pattern/glossary/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 $cref/VectorSet/RevSparseHes/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 $cref/sparsity pattern/glossary/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 $cref/elements of type/SimpleVector/Elements of Specified 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 \B{R}^{n \times n}$$ is the $latex n \times n$$ identity matrix. 00136 Further suppose that the $latex S$$ is the $th k$$ 00137 $cref/elementary vector/glossary/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 \B{R}^{n \times n}$$. 00149 00150 $head Example$$ 00151 $children% 00152 example/rev_sparse_hes.cpp 00153 %$$ 00154 The file 00155 $cref rev_sparse_hes.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 # include <cppad/local/pod_vector.hpp> 00164 # include <cppad/local/std_set.hpp> 00165 00166 CPPAD_BEGIN_NAMESPACE 00167 /*! 00168 \defgroup rev_sparse_hes_hpp rev_sparse_hes.hpp 00169 \{ 00170 \file rev_sparse_hes.hpp 00171 Reverse mode Hessian sparsity patterns. 00172 */ 00173 00174 /*! 00175 Calculate Hessian sparsity patterns using reverse mode. 00176 00177 The C++ source code corresponding to this operation is 00178 \verbatim 00179 h = f.RevSparseHes(q, s) 00180 \endverbatim 00181 00182 \tparam Base 00183 is the base type for this recording. 00184 00185 \tparam VectorSet 00186 is a simple vector with elements of type \c bool. 00187 00188 \tparam Sparsity 00189 is either \c sparse_pack, \c sparse_set, or \c sparse_list. 00190 00191 \param q 00192 is the value of \a q in the 00193 by the previous call of the form 00194 \verbatim 00195 f.ForSparseJac(q, r) 00196 \endverbatim 00197 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00198 00199 \param s 00200 is a vector with size \c m that specifies the sparsity pattern 00201 for the vector \f$ S \f$, 00202 where \c m is the number of dependent variables 00203 corresponding to the operation sequence stored in \a play. 00204 00205 \param h 00206 the input value of \a h must be a vector with size \c q*n. 00207 The input value of its elements does not matter. 00208 On output, \a h is the sparsity pattern for the matrix 00209 \f[ 00210 H(x) = R^T (S * F)^{(2)} (x) 00211 \f] 00212 where \f$ F \f$ is the function corresponding to the operation sequence 00213 and \a x is any argument value. 00214 00215 \param total_num_var 00216 is the total number of variables in this recording. 00217 00218 \param dep_taddr 00219 maps dependendent variable index 00220 to the corresponding variable in the tape. 00221 00222 \param ind_taddr 00223 maps independent variable index 00224 to the corresponding variable in the tape. 00225 00226 \param play 00227 is the recording that defines the function we are computing the sparsity 00228 pattern for. 00229 00230 \param for_jac_sparsity 00231 is a vector of sets containing the 00232 the forward Jacobian sparsity pattern corresponding to 00233 $latex R$$ for all of the variables on the tape. 00234 */ 00235 00236 template <class Base, class VectorSet, class Sparsity> 00237 void RevSparseHesBool( 00238 size_t q , 00239 const VectorSet& s , 00240 VectorSet& h , 00241 size_t total_num_var , 00242 CppAD::vector<size_t>& dep_taddr , 00243 CppAD::vector<size_t>& ind_taddr , 00244 CppAD::player<Base>& play , 00245 Sparsity& for_jac_sparsity ) 00246 { 00247 // temporary indices 00248 size_t i, j; 00249 00250 // check Vector is Simple VectorSet class with bool elements 00251 CheckSimpleVector<bool, VectorSet>(); 00252 00253 // range and domain dimensions for F 00254 size_t m = dep_taddr.size(); 00255 size_t n = ind_taddr.size(); 00256 00257 CPPAD_ASSERT_KNOWN( 00258 q == for_jac_sparsity.end(), 00259 "RevSparseHes: q (first argument) is not equal to its value\n" 00260 "in the previous call to ForSparseJac with this ADFun object." 00261 ); 00262 CPPAD_ASSERT_KNOWN( 00263 size_t(s.size()) == m, 00264 "RevSparseHes: s (second argument) length is not equal to\n" 00265 "range dimension for ADFun object." 00266 ); 00267 00268 // Array that will hold reverse Jacobian dependency flag. 00269 // Initialize as true for the dependent variables. 00270 pod_vector<bool> RevJac; 00271 RevJac.extend(total_num_var); 00272 for(i = 0; i < total_num_var; i++) 00273 RevJac[i] = false; 00274 for(i = 0; i < m; i++) 00275 { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00276 RevJac[ dep_taddr[i] ] = s[i]; 00277 } 00278 00279 00280 // vector of sets that will hold reverse Hessain values 00281 Sparsity rev_hes_sparsity; 00282 rev_hes_sparsity.resize(total_num_var, q); 00283 00284 // compute the Hessian sparsity patterns 00285 RevHesSweep( 00286 n, 00287 total_num_var, 00288 &play, 00289 for_jac_sparsity, 00290 RevJac.data(), 00291 rev_hes_sparsity 00292 ); 00293 00294 // return values corresponding to independent variables 00295 CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == n * q ); 00296 00297 // j is index corresponding to reverse mode partial 00298 for(j = 0; j < n; j++) 00299 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var ); 00300 00301 // ind_taddr[j] is operator taddr for j-th independent variable 00302 CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == j + 1 ); 00303 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00304 00305 // extract the result from rev_hes_sparsity 00306 for(i = 0; i < q; i++) 00307 h[ i * n + j ] = false; 00308 CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q ); 00309 rev_hes_sparsity.begin(j + 1); 00310 i = rev_hes_sparsity.next_element(); 00311 while( i < q ) 00312 { h[ i * n + j ] = true; 00313 i = rev_hes_sparsity.next_element(); 00314 } 00315 } 00316 00317 return; 00318 } 00319 00320 /*! 00321 Calculate Hessian sparsity patterns using reverse mode. 00322 00323 The C++ source code corresponding to this operation is 00324 \verbatim 00325 h = f.RevSparseHes(q, s) 00326 \endverbatim 00327 00328 \tparam Base 00329 is the base type for this recording. 00330 00331 \tparam VectorSet 00332 is a simple vector with elements of type \c std::set<size_t>. 00333 00334 \tparam Sparsity 00335 is either \c sparse_pack, \c sparse_set, or \c sparse_list. 00336 00337 \param q 00338 is the value of \a q in the 00339 by the previous call of the form 00340 \verbatim 00341 f.ForSparseJac(q, r) 00342 \endverbatim 00343 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00344 00345 \param s 00346 is a vector with size \c m that specifies the sparsity pattern 00347 for the vector \f$ S \f$, 00348 where \c m is the number of dependent variables 00349 corresponding to the operation sequence stored in \a play. 00350 00351 \param h 00352 the input value of \a h must be a vector with size \a q. 00353 On input, each element of \a h must be an empty set. 00354 The input value of its elements does not matter. 00355 On output, \a h is the sparsity pattern for the matrix 00356 \f[ 00357 H(x) = R^T (S * F)^{(2)} (x) 00358 \f] 00359 where \f$ F \f$ is the function corresponding to the operation sequence 00360 and \a x is any argument value. 00361 00362 \param total_num_var 00363 is the total number of variables in this recording. 00364 00365 \param dep_taddr 00366 maps dependendent variable index 00367 to the corresponding variable in the tape. 00368 00369 \param ind_taddr 00370 maps independent variable index 00371 to the corresponding variable in the tape. 00372 00373 \param play 00374 is the recording that defines the function we are computing the sparsity 00375 pattern for. 00376 00377 \param for_jac_sparsity 00378 is a vector of sets containing the 00379 the forward Jacobian sparsity pattern corresponding to 00380 $latex R$$ for all of the variables on the tape. 00381 */ 00382 00383 template <class Base, class VectorSet, class Sparsity> 00384 void RevSparseHesSet( 00385 size_t q , 00386 const VectorSet& s , 00387 VectorSet& h , 00388 size_t total_num_var , 00389 CppAD::vector<size_t>& dep_taddr , 00390 CppAD::vector<size_t>& ind_taddr , 00391 CppAD::player<Base>& play , 00392 Sparsity& for_jac_sparsity ) 00393 { 00394 // temporary indices 00395 size_t i, j; 00396 std::set<size_t>::const_iterator itr; 00397 00398 // check VectorSet is Simple Vector class with sets for elements 00399 CheckSimpleVector<std::set<size_t>, VectorSet>( 00400 one_element_std_set<size_t>(), two_element_std_set<size_t>() 00401 ); 00402 00403 // range and domain dimensions for F 00404 # ifndef NDEBUG 00405 size_t m = dep_taddr.size(); 00406 # endif 00407 size_t n = ind_taddr.size(); 00408 00409 CPPAD_ASSERT_KNOWN( 00410 q == for_jac_sparsity.end(), 00411 "RevSparseHes: q (first argument) is not equal to its value\n" 00412 "in the previous call to ForSparseJac with this ADFun object." 00413 ); 00414 CPPAD_ASSERT_KNOWN( 00415 s.size() == 1, 00416 "RevSparseHes: s (second argument) length is not equal to one." 00417 ); 00418 00419 // Array that will hold reverse Jacobian dependency flag. 00420 // Initialize as true for the dependent variables. 00421 pod_vector<bool> RevJac; 00422 RevJac.extend(total_num_var); 00423 for(i = 0; i < total_num_var; i++) 00424 RevJac[i] = false; 00425 itr = s[0].begin(); 00426 while( itr != s[0].end() ) 00427 { i = *itr++; 00428 CPPAD_ASSERT_KNOWN( 00429 i < m, 00430 "RevSparseHes: an element of the set s[0] has value " 00431 "greater than or equal m" 00432 ); 00433 CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); 00434 RevJac[ dep_taddr[i] ] = true; 00435 } 00436 00437 00438 // vector of sets that will hold reverse Hessain values 00439 Sparsity rev_hes_sparsity; 00440 rev_hes_sparsity.resize(total_num_var, q); 00441 00442 // compute the Hessian sparsity patterns 00443 RevHesSweep( 00444 n, 00445 total_num_var, 00446 &play, 00447 for_jac_sparsity, 00448 RevJac.data(), 00449 rev_hes_sparsity 00450 ); 00451 00452 // return values corresponding to independent variables 00453 // j is index corresponding to reverse mode partial 00454 CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == q ); 00455 for(j = 0; j < n; j++) 00456 { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] < total_num_var ); 00457 CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == j + 1 ); 00458 CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); 00459 00460 // extract the result from rev_hes_sparsity 00461 // and add corresponding elements to result sets in h 00462 CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q ); 00463 rev_hes_sparsity.begin(j+1); 00464 i = rev_hes_sparsity.next_element(); 00465 while( i < q ) 00466 { h[i].insert(j); 00467 i = rev_hes_sparsity.next_element(); 00468 } 00469 } 00470 00471 return; 00472 } 00473 00474 00475 00476 /*! 00477 User API for Hessian sparsity patterns using reverse mode. 00478 00479 The C++ source code corresponding to this operation is 00480 \verbatim 00481 h = f.RevSparseHes(q, r) 00482 \endverbatim 00483 00484 \tparam Base 00485 is the base type for this recording. 00486 00487 \tparam VectorSet 00488 is a simple vector with elements of type \c bool 00489 or \c std::set<size_t>. 00490 00491 \param q 00492 is the value of \a q in the 00493 by the previous call of the form 00494 \verbatim 00495 f.ForSparseJac(q, r, packed) 00496 \endverbatim 00497 The value \c r in this call is a sparsity pattern for the matrix \f$ R \f$. 00498 The type of the element of \c r for the previous call to \c ForSparseJac 00499 must be the same as the type of the elements of \c s. 00500 00501 \param s 00502 is a vector with size \c m that specifies the sparsity pattern 00503 for the vector \f$ S \f$, 00504 where \c m is the number of dependent variables 00505 corresponding to the operation sequence stored in \a play. 00506 00507 \return 00508 is a vector with size \c q*n. 00509 containing a sparsity pattern for the matrix 00510 \f[ 00511 H(x) = R^T ( S * F)^{(2)} (x) 00512 \f] 00513 where \f$ F \f$ is the function corresponding to the operation sequence 00514 and \a x is any argument value. 00515 */ 00516 00517 template <class Base> 00518 template <class VectorSet> 00519 VectorSet ADFun<Base>::RevSparseHes(size_t q, const VectorSet& s) 00520 { VectorSet h; 00521 typedef typename VectorSet::value_type Set_type; 00522 00523 // Should check to make sure q is same as in previous call to 00524 // forward sparse Jacobian. 00525 RevSparseHesCase( 00526 Set_type() , 00527 q , 00528 s , 00529 h 00530 ); 00531 00532 return h; 00533 } 00534 00535 /*! 00536 Private helper function for RevSparseHes(q, s). 00537 00538 All of the description in the public member function RevSparseHes(q, s) 00539 applies. 00540 00541 \param set_type 00542 is a \c bool value. This argument is used to dispatch to the proper source 00543 code depending on the vlaue of \c VectorSet::value_type. 00544 00545 \param q 00546 See \c RevSparseHes(q, s). 00547 00548 \param s 00549 See \c RevSparseHes(q, s). 00550 00551 \param h 00552 is the return value for the corresponging call to \c RevSparseJac(q, s). 00553 */ 00554 template <class Base> 00555 template <class VectorSet> 00556 void ADFun<Base>::RevSparseHesCase( 00557 bool set_type , 00558 size_t q , 00559 const VectorSet& s , 00560 VectorSet& h ) 00561 { size_t n = Domain(); 00562 h.resize(q * n ); 00563 00564 CPPAD_ASSERT_KNOWN( 00565 for_jac_sparse_pack_.n_set() > 0, 00566 "RevSparseHes: previous stored call to ForSparseJac did not " 00567 "use bool for the elements of r." 00568 ); 00569 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 ); 00570 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == total_num_var_ ); 00571 00572 // use sparse_pack for the calculation 00573 CppAD::RevSparseHesBool( 00574 q , 00575 s , 00576 h , 00577 total_num_var_ , 00578 dep_taddr_ , 00579 ind_taddr_ , 00580 play_ , 00581 for_jac_sparse_pack_ 00582 ); 00583 } 00584 00585 /*! 00586 Private helper function for RevSparseHes(q, s). 00587 00588 All of the description in the public member function RevSparseHes(q, s) 00589 applies. 00590 00591 \param set_type 00592 is a \c std::set<size_t> value. 00593 This argument is used to dispatch to the proper source 00594 code depending on the vlaue of \c VectorSet::value_type. 00595 00596 \param q 00597 See \c RevSparseHes(q, s). 00598 00599 \param s 00600 See \c RevSparseHes(q, s). 00601 00602 \param h 00603 is the return value for the corresponging call to \c RevSparseJac(q, s). 00604 */ 00605 template <class Base> 00606 template <class VectorSet> 00607 void ADFun<Base>::RevSparseHesCase( 00608 const std::set<size_t>& set_type , 00609 size_t q , 00610 const VectorSet& s , 00611 VectorSet& h ) 00612 { h.resize(q); 00613 00614 CPPAD_ASSERT_KNOWN( 00615 for_jac_sparse_set_.n_set() > 0, 00616 "RevSparseHes: previous stored call to ForSparseJac did not " 00617 "use std::set<size_t> for the elements of r." 00618 ); 00619 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 ); 00620 CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == total_num_var_ ); 00621 00622 // use sparse_pack for the calculation 00623 CppAD::RevSparseHesSet( 00624 q , 00625 s , 00626 h , 00627 total_num_var_ , 00628 dep_taddr_ , 00629 ind_taddr_ , 00630 play_ , 00631 for_jac_sparse_set_ 00632 ); 00633 } 00634 00635 /*! \} */ 00636 CPPAD_END_NAMESPACE 00637 # endif