CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_HESSIAN_INCLUDED 00003 # define CPPAD_SPARSE_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 sparse_hessian$$ 00018 $spell 00019 recomputed 00020 CppAD 00021 valarray 00022 std 00023 Bool 00024 hes 00025 const 00026 Taylor 00027 $$ 00028 00029 $section Sparse Hessian: Easy Driver$$ 00030 $index SparseHessian$$ 00031 $index hessian, sparse$$ 00032 00033 $head Syntax$$ 00034 $icode%hes% = %f%.SparseHessian(%x%, %w%) 00035 %hes% = %f%.SparseHessian(%x%, %w%, %p%) 00036 %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%) 00037 %$$ 00038 00039 $head Purpose$$ 00040 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size, 00041 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$. 00042 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the 00043 $cref/AD function/glossary/AD Function/$$ 00044 corresponding to $icode f$$. 00045 The syntax above sets $icode hes$$ to the Hessian 00046 $latex \[ 00047 H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 00048 \] $$ 00049 This routine takes advantage of the sparsity of the Hessian 00050 in order to reduce the amount of computation necessary. 00051 If $icode row$$ and $icode col$$ are present, it also takes 00052 advantage of the reduced set of elements of the Hessian that 00053 need to be computed. 00054 One can use speed tests (e.g. $cref speed_test$$) 00055 to verify that results are computed faster 00056 than when using the routine $cref Hessian$$. 00057 00058 $head f$$ 00059 The object $icode f$$ has prototype 00060 $codei% 00061 ADFun<%Base%> %f% 00062 %$$ 00063 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$ 00064 (see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below). 00065 00066 $head x$$ 00067 The argument $icode x$$ has prototype 00068 $codei% 00069 const %VectorBase%& %x% 00070 %$$ 00071 (see $cref/VectorBase/sparse_hessian/VectorBase/$$ below) 00072 and its size 00073 must be equal to $icode n$$, the dimension of the 00074 $cref/domain/seq_property/Domain/$$ space for $icode f$$. 00075 It specifies 00076 that point at which to evaluate the Hessian. 00077 00078 $head w$$ 00079 The argument $icode w$$ has prototype 00080 $codei% 00081 const %VectorBase%& %w% 00082 %$$ 00083 and size $latex m$$. 00084 It specifies the value of $latex w_i$$ in the expression 00085 for $icode hes$$. 00086 The more components of $latex w$$ that are identically zero, 00087 the more sparse the resulting Hessian may be (and hence the more efficient 00088 the calculation of $icode hes$$ may be). 00089 00090 $head p$$ 00091 The argument $icode p$$ is optional and has prototype 00092 $codei% 00093 const %VectorSet%& %p% 00094 %$$ 00095 (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below) 00096 If it has elements of type $code bool$$, 00097 its size is $latex n * n$$. 00098 If it has elements of type $code std::set<size_t>$$, 00099 its size is $latex n$$ and all its set elements are between 00100 zero and $latex n - 1$$. 00101 It specifies a 00102 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00103 for the Hessian $latex H(x)$$. 00104 $pre 00105 00106 $$ 00107 If this sparsity pattern does not change between calls to 00108 $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and 00109 pass this argument to $codei SparseHessian$$. 00110 In addition, 00111 if you specify $icode p$$, CppAD will use the same 00112 type of sparsity representation 00113 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00114 for its internal calculations. 00115 Otherwise, the representation 00116 for the internal calculations is unspecified. 00117 00118 $head row, col$$ 00119 The arguments $icode row$$ and $icode col$$ are optional and have prototype 00120 $codei% 00121 const %VectorSize%& %row% 00122 const %VectorSize%& %col% 00123 %$$ 00124 (see $cref/VectorSize/sparse_hessian/VectorSize/$$ below). 00125 They specify which rows and columns of $latex H (x)$$ are 00126 returned and in what order. 00127 We use $latex K$$ to denote the value $icode%hes%.size()%$$ 00128 which must also equal the size of $icode row$$ and $icode col$$. 00129 Furthermore, 00130 for $latex k = 0 , \ldots , K-1$$, it must hold that 00131 $latex row[k] < n$$ and $latex col[k] < n$$. 00132 In addition, 00133 all of the $latex (row[k], col[k])$$ pairs must correspond to a true value 00134 in the sparsity pattern $icode p$$. 00135 00136 $head hes$$ 00137 The result $icode hes$$ has prototype 00138 $codei% 00139 %VectorBase% %hes% 00140 %$$ 00141 In the case where $icode row$$ and $icode col$$ are not present, 00142 the size of $icode hes$$ is $latex n * n$$ and 00143 its size is $latex n * n$$. 00144 In this case, for $latex i = 0 , \ldots , n - 1 $$ 00145 and $latex ell = 0 , \ldots , n - 1$$ 00146 $latex \[ 00147 hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x ) 00148 \] $$ 00149 $pre 00150 00151 $$ 00152 In the case where the arguments $icode row$$ and $icode col$$ are present, 00153 we use $latex K$$ to denote the size of $icode hes$$. 00154 The input value of its elements does not matter. 00155 Upon return, for $latex k = 0 , \ldots , K - 1$$, 00156 $latex \[ 00157 hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x) 00158 \; , \; 00159 \; {\rm where} \; 00160 j = row[k] 00161 \; {\rm and } \; 00162 \ell = col[k] 00163 \] $$ 00164 00165 $head work$$ 00166 If this argument is present, it has prototype 00167 $codei% 00168 sparse_hessian_work& %work% 00169 %$$ 00170 This object can only be used with the routines $code SparseHessian$$. 00171 During its the first use, information is stored in $icode work$$. 00172 This is used to reduce the work done by future calls to $code SparseHessian$$ 00173 with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$. 00174 If a future call is make where any of these values have changed, 00175 you must first call $icode%work%.clear()%$$ 00176 to inform CppAD that this information needs to be recomputed. 00177 00178 $head n_sweep$$ 00179 The return value $icode n_sweep$$ has prototype 00180 $codei% 00181 size_t %n_sweep% 00182 %$$ 00183 It is the number of first order forward sweeps 00184 used to compute the requested Hessian values. 00185 Each first forward sweep is followed by a second order reverse sweep 00186 so it is also the number of reverse sweeps. 00187 This is proportional to the total work that $code SparseHessian$$ does, 00188 not counting the zero order forward sweep, 00189 or the work to combine multiple columns into a single 00190 forward-reverse sweep pair. 00191 00192 $head VectorBase$$ 00193 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with 00194 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00195 $icode Base$$. 00196 The routine $cref CheckSimpleVector$$ will generate an error message 00197 if this is not the case. 00198 00199 $head VectorSet$$ 00200 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with 00201 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00202 $code bool$$ or $code std::set<size_t>$$; 00203 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00204 of the difference. 00205 The routine $cref CheckSimpleVector$$ will generate an error message 00206 if this is not the case. 00207 00208 $subhead Restrictions$$ 00209 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00210 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00211 corresponding set. 00212 According to section 26.3.2.3 of the 1998 C++ standard, 00213 $code std::valarray< std::set<size_t> >$$ does not satisfy 00214 this condition. 00215 00216 $head VectorSize$$ 00217 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with 00218 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00219 $code size_t$$. 00220 The routine $cref CheckSimpleVector$$ will generate an error message 00221 if this is not the case. 00222 00223 $head Uses Forward$$ 00224 After each call to $cref Forward$$, 00225 the object $icode f$$ contains the corresponding 00226 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00227 After a call to any of the sparse Hessian routines, 00228 the zero order Taylor coefficients correspond to 00229 $icode%f%.Forward(0, %x%)%$$ 00230 and the other coefficients are unspecified. 00231 00232 $head Example$$ 00233 $children% 00234 example/sparse_hessian.cpp 00235 %$$ 00236 The routine 00237 $cref sparse_hessian.cpp$$ 00238 is examples and tests of $code sparse_hessian$$. 00239 It return $code true$$, if it succeeds and $code false$$ otherwise. 00240 00241 00242 $end 00243 ----------------------------------------------------------------------------- 00244 */ 00245 # include <cppad/local/std_set.hpp> 00246 00247 CPPAD_BEGIN_NAMESPACE 00248 /*! 00249 \defgroup sparse_hessian_hpp sparse_hessian.hpp 00250 \{ 00251 \file sparse_hessian.hpp 00252 Sparse Hessian driver routine and helper functions. 00253 */ 00254 // =========================================================================== 00255 /*! 00256 class used by SparseHessian to hold information 00257 so it does not need to be recomputed. 00258 */ 00259 class sparse_hessian_work { 00260 public: 00261 /// version of user r array sorted by row or column 00262 CppAD::vector<size_t> r_sort; 00263 /// version of user c array sorted by row or column 00264 CppAD::vector<size_t> c_sort; 00265 /// mapping from sorted array indices to user array indices 00266 CppAD::vector<size_t> k_sort; 00267 /// results of the coloring algorithm 00268 CppAD::vector<size_t> color; 00269 /// inform CppAD that this information needs to be recomputed 00270 void clear(void) 00271 { r_sort.clear(); 00272 c_sort.clear(); 00273 k_sort.clear(); 00274 color.clear(); 00275 } 00276 }; 00277 // =========================================================================== 00278 /*! 00279 Private helper function that does computation for all Sparse Hessian cases. 00280 00281 \tparam Base 00282 See \c SparseHessian(x, w, p, row, col, hes, work). 00283 00284 \tparam VectorBase 00285 See \c SparseHessian(x, w, p, row, col, hes, work). 00286 00287 \tparam VectorSet 00288 is either \c sparse_pack, \c sparse_set or \c sparse_list. 00289 00290 \param x 00291 See \c SparseHessian(x, w, p, row, col, hes, work). 00292 00293 \param w 00294 See \c SparseHessian(x, w, p, row, col, hes, work). 00295 00296 \param sparsity 00297 If <code>work.color.size() != 0</code>, then \c sparsity is not used. 00298 Otherwise, it is a 00299 sparsity pattern for the Hessian of this ADFun<Base> object. 00300 00301 \param hes 00302 See \c SparseHessian(x, w, p, row, col, hes, work). 00303 00304 \param work 00305 See \c SparseHessian(x, w, p, row, col, hes, work). 00306 00307 \return 00308 See \c SparseHessian(x, w, p, row, col, hes, work). 00309 */ 00310 template<class Base> 00311 template <class VectorBase, class VectorSet> 00312 size_t ADFun<Base>::SparseHessianCompute( 00313 const VectorBase& x , 00314 const VectorBase& w , 00315 VectorSet& sparsity , 00316 VectorBase& hes , 00317 sparse_hessian_work& work ) 00318 { 00319 using CppAD::vectorBool; 00320 size_t i, j, k, ell; 00321 00322 CppAD::vector<size_t>& row(work.r_sort); 00323 CppAD::vector<size_t>& col(work.c_sort); 00324 CppAD::vector<size_t>& user_k(work.k_sort); 00325 CppAD::vector<size_t>& color(work.color); 00326 00327 size_t n = Domain(); 00328 00329 // some values 00330 const Base zero(0); 00331 const Base one(1); 00332 00333 // check VectorBase is Simple Vector class with Base type elements 00334 CheckSimpleVector<Base, VectorBase>(); 00335 00336 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00337 CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); 00338 00339 // number of components of Hessian that are required 00340 size_t K = hes.size(); 00341 CPPAD_ASSERT_UNKNOWN( row.size() == K+1 ); 00342 CPPAD_ASSERT_UNKNOWN( col.size() == K ); 00343 CPPAD_ASSERT_UNKNOWN( row[K] == n ); 00344 00345 // Point at which we are evaluating the Hessian 00346 Forward(0, x); 00347 00348 // Rows of the Hessian (i below) correspond to the forward mode index 00349 // and columns (j below) correspond to the reverse mode index. 00350 if( color.size() == 0 ) 00351 { 00352 CPPAD_ASSERT_UNKNOWN( sparsity.n_set() == n ); 00353 CPPAD_ASSERT_UNKNOWN( sparsity.end() == n ); 00354 00355 // rows and columns that are in the returned hessian 00356 VectorSet r_used, c_used; 00357 r_used.resize(n, n); 00358 c_used.resize(n, n); 00359 k = 0; 00360 while( k < K ) 00361 { CPPAD_ASSERT_UNKNOWN( row[k] < n && col[k] < n ); 00362 CPPAD_ASSERT_UNKNOWN( k == 0 || row[k-1] <= row[k] ); 00363 CPPAD_ASSERT_KNOWN( 00364 sparsity.is_element(row[k], col[k]) , 00365 "SparseHessian: an (row, col) pair is not in sparsity pattern." 00366 ); 00367 r_used.add_element(col[k], row[k]); 00368 c_used.add_element(row[k], col[k]); 00369 k++; 00370 } 00371 00372 // given a column index, which rows are non-zero and not used 00373 VectorSet not_used; 00374 not_used.resize(n, n); 00375 for(i = 0; i < n; i++) 00376 { sparsity.begin(i); 00377 j = sparsity.next_element(); 00378 while( j != sparsity.end() ) 00379 { if( ! r_used.is_element(j , i) ) 00380 not_used.add_element(j, i); 00381 j = sparsity.next_element(); 00382 } 00383 } 00384 00385 // initial coloring 00386 color.resize(n); 00387 for(i = 0; i < n; i++) 00388 color[i] = i; 00389 00390 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00391 // Graph Coloring in Optimization Revisited by 00392 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00393 vectorBool forbidden(n); 00394 for(i = 1; i < n; i++) 00395 { 00396 // initial all colors as ok for this row 00397 // (value of forbidden for ell > i does not matter) 00398 for(ell = 0; ell <= i; ell++) 00399 forbidden[ell] = false; 00400 00401 // ----------------------------------------------------- 00402 // Forbid colors that this row would destroy results for. 00403 // for each column that is non-zero for this row 00404 sparsity.begin(i); 00405 j = sparsity.next_element(); 00406 while( j != sparsity.end() ) 00407 { // for each row that this column uses 00408 r_used.begin(j); 00409 ell = r_used.next_element(); 00410 while( ell != r_used.end() ) 00411 { // if this is not the same row, forbid its color 00412 if( ell < i ) 00413 forbidden[ color[ell] ] = true; 00414 ell = r_used.next_element(); 00415 } 00416 j = sparsity.next_element(); 00417 } 00418 00419 // ------------------------------------------------------- 00420 // Forbid colors that would destroy the results for this row. 00421 // for each column that this row used 00422 c_used.begin(i); 00423 j = c_used.next_element(); 00424 while( j != c_used.end() ) 00425 { // For each row that is non-zero for this column 00426 // (the used rows have already been checked above). 00427 not_used.begin(j); 00428 ell = not_used.next_element(); 00429 while( ell != not_used.end() ) 00430 { // if this is not the same row, forbid its color 00431 if( ell < i ) 00432 forbidden[ color[ell] ] = true; 00433 ell = not_used.next_element(); 00434 } 00435 j = c_used.next_element(); 00436 } 00437 00438 // pick the color with the smallest index 00439 ell = 0; 00440 while( forbidden[ell] ) 00441 { ell++; 00442 CPPAD_ASSERT_UNKNOWN( ell <= i ); 00443 } 00444 color[i] = ell; 00445 } 00446 } 00447 size_t n_color = 1; 00448 for(ell = 0; ell < n; ell++) 00449 n_color = std::max(n_color, color[ell] + 1); 00450 00451 // direction vector for calls to forward (rows of the Hessian) 00452 VectorBase u(n); 00453 00454 // location for return values from reverse (columns of the Hessian) 00455 VectorBase ddw(2 * n); 00456 00457 // initialize the return value 00458 for(k = 0; k < K; k++) 00459 hes[k] = zero; 00460 00461 // loop over colors 00462 size_t n_sweep = 0; 00463 for(ell = 0; ell < n_color; ell++) 00464 { bool any = false; 00465 k = 0; 00466 for(i = 0; i < n; i++) if( color[i] == ell ) 00467 { // find first k such that row[k] has color ell 00468 if( ! any ) 00469 { while( row[k] < i ) 00470 k++; 00471 any = row[k] == i; 00472 } 00473 } 00474 if( any ) 00475 { n_sweep++; 00476 // combine all rows with this color 00477 for(i = 0; i < n; i++) 00478 { u[i] = zero; 00479 if( color[i] == ell ) 00480 u[i] = one; 00481 } 00482 // call forward mode for all these rows at once 00483 Forward(1, u); 00484 00485 // evaluate derivative of w^T * F'(x) * u 00486 ddw = Reverse(2, w); 00487 00488 // set the corresponding components of the result 00489 for(i = 0; i < n; i++) if( color[i] == ell ) 00490 { // find first index in c for this column 00491 while( row[k] < i ) 00492 k++; 00493 // extract the results for this row 00494 while( row[k] == i ) 00495 { hes[ user_k[k] ] = ddw[ col[k] * 2 + 1 ]; 00496 k++; 00497 } 00498 } 00499 } 00500 } 00501 return n_sweep; 00502 } 00503 // =========================================================================== 00504 /*! 00505 Private helper function for vector of \c bool sparsity pattern cases. 00506 00507 \tparam Base 00508 See \c SparseHessian(x, w, p, row, col, hes, work). 00509 00510 \tparam VectorBase 00511 See \c SparseHessian(x, w, p, row, col, hes, work). 00512 00513 \tparam VectorSet 00514 is a simple vector with elements of type \c bool. 00515 00516 \param set_type 00517 has element type for vector representing the sparsity sets. 00518 00519 \param x 00520 See \c SparseHessian(x, w, p, row, col, hes, work). 00521 00522 \param w 00523 See \c SparseHessian(x, w, p, row, col, hes, work). 00524 00525 \param p 00526 Sparsity pattern for the Hessian of this ADFun<Base> object. 00527 00528 \param hes 00529 See \c SparseHessian(x, w, p, row, col, hes, work). 00530 00531 \param work 00532 See \c SparseHessian(x, w, p, row, col, hes, work). 00533 00534 \return 00535 See \c SparseHessian(x, w, p, row, col, hes, work). 00536 */ 00537 template <class Base> 00538 template <class VectorBase, class VectorSet> 00539 size_t ADFun<Base>::SparseHessianCase( 00540 bool set_type , 00541 const VectorBase& x , 00542 const VectorBase& w , 00543 const VectorSet& p , 00544 VectorBase& hes , 00545 sparse_hessian_work& work ) 00546 { 00547 size_t n = Domain(); 00548 00549 // check VectorSet is Simple Vector class with bool elements 00550 CheckSimpleVector<bool, VectorSet>(); 00551 00552 // check VectorBase is Simple Vector class with Base type elements 00553 CheckSimpleVector<Base, VectorBase>(); 00554 00555 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00556 CPPAD_ASSERT_UNKNOWN( size_t(w.size()) == Range() ); 00557 CPPAD_ASSERT_KNOWN( 00558 size_t(p.size()) == n * n, 00559 "SparseHessian: using bool values for sparsity and p.size() " 00560 "not equal square of domain dimension for f" 00561 ); 00562 00563 sparse_pack sparsity; 00564 if( work.color.size() == 0 ) 00565 { bool transpose = false; 00566 typedef typename VectorSet::value_type Set_type; 00567 sparsity_user2internal(sparsity, p, n, n, transpose); 00568 } 00569 00570 // compute the Hessian 00571 size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work); 00572 00573 return n_sweep; 00574 } 00575 /*! 00576 Private helper function for vector of std::set<size_t> sparsity pattern cases. 00577 00578 \tparam Base 00579 See \c SparseHessian(x, w, p, row, col, hes, work). 00580 00581 \tparam VectorBase 00582 See \c SparseHessian(x, w, p, row, col, hes, work). 00583 00584 \tparam VectorSet 00585 is a simple vector with elements of type <code>std::set<size_t></code>. 00586 00587 \param set_type 00588 has element type for vector representing the sparsity sets. 00589 00590 \param x 00591 See \c SparseHessian(x, w, p, row, col, hes, work). 00592 00593 \param w 00594 See \c SparseHessian(x, w, p, row, col, hes, work). 00595 00596 \param p 00597 Sparsity pattern for the Hessian of this ADFun<Base> object. 00598 00599 \param hes 00600 See \c SparseHessian(x, w, p, row, col, hes, work). 00601 00602 \param work 00603 See \c SparseHessian(x, w, p, row, col, hes, work). 00604 00605 \return 00606 See \c SparseHessian(x, w, p, row, col, hes, work). 00607 */ 00608 template <class Base> 00609 template <class VectorBase, class VectorSet> 00610 size_t ADFun<Base>::SparseHessianCase( 00611 const std::set<size_t>& set_type , 00612 const VectorBase& x , 00613 const VectorBase& w , 00614 const VectorSet& p , 00615 VectorBase& hes , 00616 sparse_hessian_work& work ) 00617 { 00618 size_t n = Domain(); 00619 00620 // check VectorSet is Simple Vector class with std::set<size_t> elements 00621 CheckSimpleVector<std::set<size_t>, VectorSet>( 00622 one_element_std_set<size_t>(), two_element_std_set<size_t>() 00623 ); 00624 00625 // check VectorBase is Simple Vector class with Base type elements 00626 CheckSimpleVector<Base, VectorBase>(); 00627 00628 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00629 CPPAD_ASSERT_UNKNOWN( size_t(w.size()) == Range() ); 00630 CPPAD_ASSERT_KNOWN( 00631 size_t(p.size()) == n, 00632 "SparseHessian: using std::set<size_t> for sparsity and p.size() " 00633 "not equal domain dimension for f" 00634 ); 00635 00636 CPPAD_INTERNAL_SPARSE_SET sparsity; 00637 if( work.color.size() == 0 ) 00638 { bool transpose = false; 00639 typedef typename VectorSet::value_type Set_type; 00640 sparsity_user2internal(sparsity, p, n, n, transpose); 00641 } 00642 00643 // compute the Hessian 00644 size_t n_sweep = SparseHessianCompute(x, w, sparsity, hes, work); 00645 00646 return n_sweep; 00647 } 00648 // =========================================================================== 00649 /*! 00650 Private helper function for SparseHessian(x, w, p). 00651 00652 All of the description in the public member function SparseHessian(x, w, p) 00653 applies. 00654 00655 \param set_type 00656 is a \c bool value. This argument is used to dispatch to the proper source 00657 code depending on the value of \c VectorSet::value_type. 00658 00659 \param x 00660 See \c SparseHessian(x, w, p). 00661 00662 \param w 00663 See \c SparseHessian(x, w, p). 00664 00665 \param p 00666 See \c SparseHessian(x, w, p). 00667 00668 \param hes 00669 is the return value for the corresponding call to \c SparseHessian(x, w, p). 00670 On input, it must have size equal to the domain times range dimension 00671 for this ADFun<Base> object. 00672 On return, it will contain the Hessian. 00673 */ 00674 00675 template <class Base> 00676 template <class VectorBase, class VectorSet> 00677 void ADFun<Base>::SparseHessianCase( 00678 bool set_type , 00679 const VectorBase& x , 00680 const VectorBase& w , 00681 const VectorSet& p , 00682 VectorBase& hes ) 00683 { 00684 size_t n = Domain(); 00685 size_t i, j, k; 00686 00687 // check Vector is Simple Vector class with bool elements 00688 CheckSimpleVector<bool, VectorSet>(); 00689 00690 // check Vector is Simple Vector class with Base type elements 00691 CheckSimpleVector<Base, VectorBase>(); 00692 00693 CPPAD_ASSERT_KNOWN( 00694 size_t(x.size()) == n, 00695 "SparseHessian: size of x not equal domain dimension for f" 00696 ); 00697 CPPAD_ASSERT_KNOWN( 00698 size_t(w.size()) == Range(), 00699 "SparseHessian: size of w not equal range dimension for f" 00700 ); 00701 CPPAD_ASSERT_KNOWN( 00702 size_t(p.size()) == n * n, 00703 "SparseHessian: using bool values and size of p " 00704 "not equal square of domain dimension for f" 00705 ); 00706 sparse_hessian_work work; 00707 CppAD::vector<size_t>& row(work.r_sort); 00708 CppAD::vector<size_t>& col(work.c_sort); 00709 CppAD::vector<size_t>& user_k(work.k_sort); 00710 00711 k = 0; 00712 for(i = 0; i < n; i++) 00713 { for(j = 0; j < n; j++) 00714 if( p[i * n + j] ) 00715 k++; 00716 } 00717 size_t K = k; 00718 VectorBase H(K); 00719 row.resize(K+1); 00720 col.resize(K); 00721 user_k.resize(K); 00722 k = 0; 00723 for(i = 0; i < n; i++) 00724 { for(j = 0; j < n; j++) 00725 if( p[i * n + j] ) 00726 { row[k] = i; 00727 col[k] = j; 00728 user_k[k] = k; 00729 k++; 00730 } 00731 } 00732 row[K] = n; 00733 00734 SparseHessianCase(set_type, x, w, p, H, work); 00735 00736 Base zero(0); 00737 for(i = 0; i < n; i++) 00738 { for(j = 0; j < n; j++) 00739 hes[i * n + j] = zero; 00740 } 00741 for(k = 0; k < K; k++) 00742 hes[row[k] * n + col[k]] = H[k]; 00743 } 00744 /*! 00745 Private helper function for SparseHessian(x, w, p). 00746 00747 All of the description in the public member function SparseHessian(x, w, p) 00748 applies. 00749 00750 \param set_type 00751 is a \c std::set<size_t> value. 00752 This argument is used to dispatch to the proper source 00753 code depending on the value of \c VectorSet::value_type. 00754 00755 \param x 00756 See \c SparseHessian(x, w, p). 00757 00758 \param w 00759 See \c SparseHessian(x, w, p). 00760 00761 \param p 00762 See \c SparseHessian(x, w, p). 00763 00764 \param hes 00765 is the return value for the corresponding call to \c SparseHessian(x, w, p). 00766 On input, it must have size equal to the domain times range dimension 00767 for this ADFun<Base> object. 00768 On return, it will contain the Hessian. 00769 */ 00770 00771 template <class Base> 00772 template <class VectorBase, class VectorSet> 00773 void ADFun<Base>::SparseHessianCase( 00774 const std::set<size_t>& set_type , 00775 const VectorBase& x , 00776 const VectorBase& w , 00777 const VectorSet& p , 00778 VectorBase& hes ) 00779 { 00780 size_t n = Domain(); 00781 size_t i, j, k; 00782 00783 // check VectorSet is Simple Vector class with sets for elements 00784 CheckSimpleVector<std::set<size_t>, VectorSet>( 00785 one_element_std_set<size_t>(), two_element_std_set<size_t>() 00786 ); 00787 00788 // check Vector is Simple Vector class with Base type elements 00789 CheckSimpleVector<Base, VectorBase>(); 00790 00791 CPPAD_ASSERT_KNOWN( 00792 size_t(x.size()) == n, 00793 "SparseHessian: size of x not equal domain dimension for f" 00794 ); 00795 CPPAD_ASSERT_KNOWN( 00796 size_t(w.size()) == Range(), 00797 "SparseHessian: size of w not equal range dimension for f" 00798 ); 00799 CPPAD_ASSERT_KNOWN( 00800 size_t(p.size()) == n, 00801 "SparseHessian: using size_t sets for sparsity pattern and p.size() " 00802 "not equal domain dimension for f" 00803 ); 00804 sparse_hessian_work work; 00805 CppAD::vector<size_t>& row(work.r_sort); 00806 CppAD::vector<size_t>& col(work.c_sort); 00807 CppAD::vector<size_t>& user_k(work.k_sort); 00808 00809 k = 0; 00810 std::set<size_t>::const_iterator itr; 00811 for(i = 0; i < n; i++) 00812 { itr = p[i].begin(); 00813 while( itr != p[i].end() ) 00814 { itr++; 00815 k++; 00816 } 00817 } 00818 size_t K = k; 00819 VectorBase H(K); 00820 row.resize(K+1); 00821 col.resize(K); 00822 user_k.resize(K); 00823 k = 0; 00824 for(i = 0; i < n; i++) 00825 { itr = p[i].begin(); 00826 while( itr != p[i].end() ) 00827 { row[k] = i; 00828 col[k] = *itr; 00829 user_k[k] = k; 00830 itr++; 00831 k++; 00832 } 00833 } 00834 row[K] = n; 00835 00836 SparseHessianCase(set_type, x, w, p, H, work); 00837 00838 Base zero(0); 00839 for(i = 0; i < n; i++) 00840 { for(j = 0; j < n; j++) 00841 hes[i * n + j] = zero; 00842 } 00843 for(k = 0; k < K; k++) 00844 hes[row[k] * n + col[k]] = H[k]; 00845 } 00846 // =========================================================================== 00847 // Public Member Functions 00848 // =========================================================================== 00849 /*! 00850 Compute user specified subset of a sparse Hessian. 00851 00852 The C++ source code corresponding to this operation is 00853 \verbatim 00854 SparceHessian(x, w, p, row, col, hes, work) 00855 \endverbatim 00856 00857 \tparam Base 00858 is the base type for the recording that is stored in this ADFun<Base object. 00859 00860 \tparam VectorBase 00861 is a simple vector class with elements of type \a Base. 00862 00863 \tparam VectorSet 00864 is a simple vector class with elements of type 00865 \c bool or \c std::set<size_t>. 00866 00867 \tparam VectorSize 00868 is a simple vector class with elements of type \c size_t. 00869 00870 \param x 00871 is a vector specifing the point at which to compute the Hessian. 00872 00873 \param w 00874 is the weighting vector that defines a scalar valued function by 00875 a weighted sum of the components of the vector valued function 00876 $latex F(x)$$. 00877 00878 \param p 00879 is the sparsity pattern for the Hessian that we are calculating. 00880 00881 \param row 00882 is the vector of row indices for the returned Hessian values. 00883 00884 \param col 00885 is the vector of columns indices for the returned Hessian values. 00886 It must have the same size are r. 00887 00888 \param hes 00889 is the vector of Hessian values. 00890 It must have the same size are r. 00891 The return value <code>hes[k]</code> is the second partial of 00892 \f$ w^{\rm T} F(x)\f$ with respect to the 00893 <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. 00894 00895 \param work 00896 contains information that depends on the function object, sparsity pattern, 00897 \c row, and \c col vector. 00898 If these values are the same, \c work does not need to be recomputed. 00899 To be more specific, 00900 \c r_sort is sorted copy of \c row , 00901 \c c_sort is sorted copy of \c col , 00902 <code>k_sort[k]</code> is the original index corresponding to the 00903 values <code>r_sort[k]</code> and <code>c_sort[k]</code>. 00904 The order for the sort is by columns. 00905 Let \c n the domain dimension, 00906 and \c K the size of \c row , \c col , and \c hes. 00907 There is one extra entry 00908 in the sorted row array and it has value <code>r_sort[K]=n</code>. 00909 The \c color vector is set and used by \c SparseHessianCompute. 00910 00911 \return 00912 Is the number of first order forward sweeps used to compute the 00913 requested Hessian values. 00914 (This is also equal to the number of second order reverse sweeps.) 00915 The total work, not counting the zero order 00916 forward sweep, or the time to combine computations, is proportional to this 00917 return value. 00918 */ 00919 template<class Base> 00920 template <class VectorBase, class VectorSet, class VectorSize> 00921 size_t ADFun<Base>::SparseHessian( 00922 const VectorBase& x , 00923 const VectorBase& w , 00924 const VectorSet& p , 00925 const VectorSize& row , 00926 const VectorSize& col , 00927 VectorBase& hes , 00928 sparse_hessian_work& work ) 00929 { 00930 size_t n = Domain(); 00931 size_t k, K = hes.size(); 00932 if( work.r_sort.size() == 0 ) 00933 { // create version of (row, col, k) sorted by row value 00934 work.c_sort.resize(K); 00935 work.r_sort.resize(K+1); 00936 work.k_sort.resize(K); 00937 00938 // put sorting indices in k_sort 00939 index_sort(row, work.k_sort); 00940 00941 for(k = 0; k < K; k++) 00942 { work.r_sort[k] = row[ work.k_sort[k] ]; 00943 work.c_sort[k] = col[ work.k_sort[k] ]; 00944 } 00945 work.r_sort[K] = n; 00946 } 00947 # ifndef NDEBUG 00948 CPPAD_ASSERT_KNOWN( 00949 size_t(x.size()) == n , 00950 "SparseHessian: size of x not equal domain dimension for f." 00951 ); 00952 CPPAD_ASSERT_KNOWN( 00953 size_t(row.size()) == K && size_t(col.size()) == K , 00954 "SparseHessian: either r or c does not have the same size as ehs." 00955 ); 00956 CPPAD_ASSERT_KNOWN( 00957 work.r_sort.size() == K+1 && 00958 work.c_sort.size() == K && 00959 work.k_sort.size() == K , 00960 "SparseHessian: invalid value in work." 00961 ); 00962 CPPAD_ASSERT_KNOWN( 00963 work.color.size() == 0 || work.color.size() == n, 00964 "SparseHessian: invalid value in work." 00965 ); 00966 for(k = 0; k < K; k++) 00967 { CPPAD_ASSERT_KNOWN( 00968 row[k] < n, 00969 "SparseHessian: invalid value in r." 00970 ); 00971 CPPAD_ASSERT_KNOWN( 00972 col[k] < n, 00973 "SparseHessian: invalid value in c." 00974 ); 00975 CPPAD_ASSERT_KNOWN( 00976 work.k_sort[k] < K, 00977 "SparseHessian: invalid value in work." 00978 ); 00979 CPPAD_ASSERT_KNOWN( 00980 work.r_sort[k] == row[ work.k_sort[k] ] , 00981 "SparseHessian: invalid value in work." 00982 ); 00983 CPPAD_ASSERT_KNOWN( 00984 work.c_sort[k] == col[ work.k_sort[k] ], 00985 "SparseHessian: invalid value in work." 00986 ); 00987 } 00988 if( work.color.size() != 0 ) 00989 for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN( 00990 work.color[j] < n, 00991 "SparseHessian: invalid value in work." 00992 ); 00993 # endif 00994 00995 typedef typename VectorSet::value_type Set_type; 00996 size_t n_sweep = SparseHessianCase(Set_type(), x, w, p, hes, work); 00997 return n_sweep; 00998 } 00999 /*! 01000 Compute a sparse Hessian. 01001 01002 The C++ source code coresponding to this operation is 01003 \verbatim 01004 hes = SparseHessian(x, w, p) 01005 \endverbatim 01006 01007 01008 \tparam Base 01009 is the base type for the recording that is stored in this 01010 ADFun<Base object. 01011 01012 \tparam VectorBase 01013 is a simple vector class with elements of the \a Base. 01014 01015 \tparam VectorSet 01016 is a simple vector class with elements of type 01017 \c bool or \c std::set<size_t>. 01018 01019 \param x 01020 is a vector specifing the point at which to compute the Hessian. 01021 01022 \param w 01023 The Hessian is computed for a weighted sum of the components 01024 of the function corresponding to this ADFun<Base> object. 01025 The argument \a w specifies the weights for each component. 01026 It must have size equal to the range dimension for this ADFun<Base> object. 01027 01028 \param p 01029 is a sparsity pattern for the Hessian. 01030 01031 \return 01032 Will be a vector of size \c n * n containing the Hessian of 01033 at the point specified by \a x 01034 (where \c n is the domain dimension for this ADFun<Base> object). 01035 */ 01036 template <class Base> 01037 template <class VectorBase, class VectorSet> 01038 VectorBase ADFun<Base>::SparseHessian( 01039 const VectorBase& x, const VectorBase& w, const VectorSet& p 01040 ) 01041 { size_t n = Domain(); 01042 VectorBase hes(n * n); 01043 01044 typedef typename VectorSet::value_type Set_type; 01045 SparseHessianCase(Set_type(), x, w, p, hes); 01046 01047 return hes; 01048 } 01049 /*! 01050 Compute a sparse Hessian 01051 01052 The C++ source code coresponding to this operation is 01053 \verbatim 01054 hes = SparseHessian(x, w) 01055 \endverbatim 01056 01057 01058 \tparam Base 01059 is the base type for the recording that is stored in this 01060 ADFun<Base object. 01061 01062 \tparam VectorBase 01063 is a simple vector class with elements of the \a Base. 01064 01065 \param x 01066 is a vector specifing the point at which to compute the Hessian. 01067 01068 \param w 01069 The Hessian is computed for a weighted sum of the components 01070 of the function corresponding to this ADFun<Base> object. 01071 The argument \a w specifies the weights for each component. 01072 It must have size equal to the range dimension for this ADFun<Base> object. 01073 01074 \return 01075 Will be a vector of size \c n * n containing the Hessian of 01076 at the point specified by \a x 01077 (where \c n is the domain dimension for this ADFun<Base> object). 01078 */ 01079 template <class Base> 01080 template <class VectorBase> 01081 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w) 01082 { size_t i, j, k; 01083 typedef CppAD::vectorBool VectorBool; 01084 01085 size_t m = Range(); 01086 size_t n = Domain(); 01087 01088 // determine the sparsity pattern p for Hessian of w^T F 01089 VectorBool r(n * n); 01090 for(j = 0; j < n; j++) 01091 { for(k = 0; k < n; k++) 01092 r[j * n + k] = false; 01093 r[j * n + j] = true; 01094 } 01095 ForSparseJac(n, r); 01096 // 01097 VectorBool s(m); 01098 for(i = 0; i < m; i++) 01099 s[i] = w[i] != 0; 01100 VectorBool p = RevSparseHes(n, s); 01101 01102 // compute sparse Hessian 01103 VectorBase hes(n * n); 01104 bool set_type = true; // only used to dispatch compiler to proper case 01105 SparseHessianCase(set_type, x, w, p, hes); 01106 01107 return hes; 01108 } 01109 01110 /*! \} */ 01111 CPPAD_END_NAMESPACE 01112 # endif