CppAD: A C++ Algorithmic Differentiation Package  20171217
sparse_hessian.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin sparse_hessian$$17 spell 18 jacobian 19 recomputed 20 CppAD 21 valarray 22 std 23 Bool 24 hes 25 const 26 Taylor 27 cppad 28 cmake 29 colpack 30$$ 31 32$section Sparse Hessian$$33 34 head Syntax$$
35 $icode%hes% = %f%.SparseHessian(%x%, %w%) 36 %hes% = %f%.SparseHessian(%x%, %w%, %p%) 37 %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%) 38 %$$39 40 head Purpose$$ 41 We use$latex n$$for the cref/domain/seq_property/Domain/$$ size,
42 and $latex m$$for the cref/range/seq_property/Range/$$ size of$icode f$$. 43 We use latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
44 $cref/AD function/glossary/AD Function/$$45 corresponding to icode f$$. 46 The syntax above sets$icode hes$$to the Hessian 47 latex $48 H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 49$$$
50 This routine takes advantage of the sparsity of the Hessian
51 in order to reduce the amount of computation necessary.
52 If $icode row$$and icode col$$ are present, it also takes 53 advantage of the reduced set of elements of the Hessian that 54 need to be computed. 55 One can use speed tests (e.g.$cref speed_test$$) 56 to verify that results are computed faster 57 than when using the routine cref Hessian$$.
58
59 $head f$$60 The object icode f$$ has prototype 61$codei%
63 %$$64 Note that the cref ADFun$$ object $icode f$$is not code const$$ 65 (see$cref/Uses Forward/sparse_hessian/Uses Forward/$$below). 66 67 head x$$
68 The argument $icode x$$has prototype 69 codei% 70 const %VectorBase%& %x% 71 %$$ 72 (see$cref/VectorBase/sparse_hessian/VectorBase/$$below) 73 and its size 74 must be equal to icode n$$, the dimension of the
75 $cref/domain/seq_property/Domain/$$space for icode f$$. 76 It specifies 77 that point at which to evaluate the Hessian. 78 79$head w$$80 The argument icode w$$ has prototype
81 $codei% 82 const %VectorBase%& %w% 83 %$$84 and size latex m$$. 85 It specifies the value of$latex w_i$$in the expression 86 for icode hes$$.
87 The more components of $latex w$$that are identically zero, 88 the more sparse the resulting Hessian may be (and hence the more efficient 89 the calculation of icode hes$$ may be). 90 91$head p$$92 The argument icode p$$ is optional and has prototype
93 $codei% 94 const %VectorSet%& %p% 95 %$$96 (see cref/VectorSet/sparse_hessian/VectorSet/$$ below) 97 If it has elements of type$code bool$$, 98 its size is latex n * n$$.
99 If it has elements of type $code std::set<size_t>$$, 100 its size is latex n$$ and all its set elements are between 101 zero and$latex n - 1$$. 102 It specifies a 103 cref/sparsity pattern/glossary/Sparsity Pattern/$$
104 for the Hessian $latex H(x)$$. 105 106 subhead Purpose$$ 107 If this sparsity pattern does not change between calls to 108$codei SparseHessian$$, it should be faster to calculate icode p$$ once and
109 pass this argument to $codei SparseHessian$$. 110 If you specify icode p$$, CppAD will use the same 111 type of sparsity representation 112 (vectors of$code bool$$or vectors of code std::set<size_t>$$)
113 for its internal calculations.
114 Otherwise, the representation
115 for the internal calculations is unspecified.
116
117 $subhead work$$118 If you specify icode work$$ in the calling sequence, 119 it is not necessary to keep the sparsity pattern; see the heading 120$cref/p/sparse_hessian/work/p/$$under the icode work$$ description.
121
122 $subhead Column Subset$$123 If the arguments icode row$$ and$icode col$$are present, 124 and cref/color_method/sparse_hessian/work/color_method/$$ is
125 $code cppad.general$$or code cppad.symmetric$$, 126 it is not necessary to compute the entire sparsity pattern. 127 Only the following subset of column values will matter: 128$codei%
129  { %col%[%k%] : %k% = 0 , %...% , %K%-1 }
130 %$$. 131 132 133 head row, col$$
134 The arguments $icode row$$and icode col$$ are optional and have prototype 135$codei%
136  const %VectorSize%& %row%
137  const %VectorSize%& %col%
138 %$$139 (see cref/VectorSize/sparse_hessian/VectorSize/$$ below).
140 They specify which rows and columns of $latex H (x)$$are 141 returned and in what order. 142 We use latex K$$ to denote the value$icode%hes%.size()%$$143 which must also equal the size of icode row$$ and $icode col$$. 144 Furthermore, 145 for latex k = 0 , \ldots , K-1$$, it must hold that 146$latex row[k] < n$$and latex col[k] < n$$.
148 all of the $latex (row[k], col[k])$$pairs must correspond to a true value 149 in the sparsity pattern icode p$$. 150 151$head hes$$152 The result icode hes$$ has prototype
153 $codei% 154 %VectorBase% %hes% 155 %$$156 In the case where icode row$$ and$icode col$$are not present, 157 the size of icode hes$$ is $latex n * n$$and 158 its size is latex n * n$$. 159 In this case, for$latex i = 0 , \ldots , n - 1 $$160 and latex ell = 0 , \ldots , n - 1$$
161 $latex $162 hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x ) 163$ $$164 pre 165 166$$ 167 In the case where the arguments$icode row$$and icode col$$ are present,
168 we use $latex K$$to denote the size of icode hes$$. 169 The input value of its elements does not matter. 170 Upon return, for$latex k = 0 , \ldots , K - 1$$, 171 latex $172 hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x) 173 \; , \; 174 \; {\rm where} \; 175 j = row[k] 176 \; {\rm and } \; 177 \ell = col[k] 178$$$
179
180 $head work$$181 If this argument is present, it has prototype 182 codei% 183 sparse_hessian_work& %work% 184 %$$ 185 This object can only be used with the routines$code SparseHessian$$. 186 During its the first use, information is stored in icode work$$.
187 This is used to reduce the work done by future calls to $code SparseHessian$$188 with the same icode f$$,$icode p$$, icode row$$, and $icode col$$. 189 If a future call is made where any of these values have changed, 190 you must first call icode%work%.clear()%$$ 191 to inform CppAD that this information needs to be recomputed. 192 193$subhead color_method$$194 The coloring algorithm determines which rows and columns 195 can be computed during the same sweep. 196 This field has prototype 197 codei% 198 std::string %work%.color_method 199 %$$
200 This value only matters on the first call to $code sparse_hessian$$that 201 follows the icode work$$ constructor or a call to 202$icode%work%.clear()%$$. 203 codei% 204 205 "cppad.symmetric" 206 %$$
207 This is the default coloring method (after a constructor or $code clear()$$). 208 It takes advantage of the fact that the Hessian matrix 209 is symmetric to find a coloring that requires fewer 210 cref/sweeps/sparse_hessian/n_sweep/$$. 211$codei%
212
214 %$$215 This is the same as the code "cppad"$$ method for the
216 $cref/sparse_jacobian/sparse_jacobian/work/color_method/$$calculation. 217 codei% 218 219 "colpack.symmetric" 220 %$$ 221 This method requires that 222$cref colpack_prefix$$was specified on the 223 cref/cmake command/cmake/CMake Command/$$ line.
224 It also takes advantage of the fact that the Hessian matrix is symmetric.
225 $codei% 226 227 "colpack.general" 228 %$$229 This is the same as the code "colpack"$$ method for the 230$cref/sparse_jacobian/sparse_jacobian/work/color_method/$$calculation. 231 232 subhead colpack.star Deprecated 2017-06-01$$
233 The $code colpack.star$$method is deprecated. 234 It is the same as the code colpack.symmetric$$ 235 which should be used instead. 236 237$subhead p$$238 If icode work$$ is present, and it is not the first call after
239 its construction or a clear,
240 the sparsity pattern $icode p$$is not used. 241 This enables one to free the sparsity pattern 242 and still compute corresponding sparse Hessians. 243 244 head n_sweep$$ 245 The return value$icode n_sweep$$has prototype 246 codei% 247 size_t %n_sweep% 248 %$$
249 It is the number of first order forward sweeps
250 used to compute the requested Hessian values.
251 Each first forward sweep is followed by a second order reverse sweep
252 so it is also the number of reverse sweeps.
253 This is proportional to the total work that $code SparseHessian$$does, 254 not counting the zero order forward sweep, 255 or the work to combine multiple columns into a single 256 forward-reverse sweep pair. 257 258 head VectorBase$$ 259 The type$icode VectorBase$$must be a cref SimpleVector$$ class with
260 $cref/elements of type/SimpleVector/Elements of Specified Type/$$261 icode Base$$. 262 The routine$cref CheckSimpleVector$$will generate an error message 263 if this is not the case. 264 265 head VectorSet$$
266 The type $icode VectorSet$$must be a cref SimpleVector$$ class with 267$cref/elements of type/SimpleVector/Elements of Specified Type/$$268 code bool$$ or $code std::set<size_t>$$; 269 see cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 270 of the difference. 271 The routine$cref CheckSimpleVector$$will generate an error message 272 if this is not the case. 273 274 subhead Restrictions$$
275 If $icode VectorSet$$has elements of code std::set<size_t>$$, 276 then$icode%p%[%i%]%$$must return a reference (not a copy) to the 277 corresponding set. 278 According to section 26.3.2.3 of the 1998 C++ standard, 279 code std::valarray< std::set<size_t> >$$ does not satisfy
280 this condition.
281
282 $head VectorSize$$283 The type icode VectorSize$$ must be a$cref SimpleVector$$class with 284 cref/elements of type/SimpleVector/Elements of Specified Type/$$
285 $code size_t$$. 286 The routine cref CheckSimpleVector$$ will generate an error message 287 if this is not the case. 288 289$head Uses Forward$$290 After each call to cref Forward$$,
291 the object $icode f$$contains the corresponding 292 cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 293 After a call to any of the sparse Hessian routines, 294 the zero order Taylor coefficients correspond to 295$icode%f%.Forward(0, %x%)%$$296 and the other coefficients are unspecified. 297 298 children% 299 example/sparse/sparse_hessian.cpp% 300 example/sparse/sub_sparse_hes.cpp% 301 example/sparse/sparse_sub_hes.cpp 302 %$$
303
304 $head Example$$305 The routine 306 cref sparse_hessian.cpp$$ 307 is examples and tests of$code sparse_hessian$$. 308 It return code true$$, if it succeeds and $code false$$otherwise. 309 310 head Subset Hessian$$ 311 The routine 312$cref sub_sparse_hes.cpp$$313 is an example and test that compute a sparse Hessian 314 for a subset of the variables. 315 It returns code true$$, for success, and $code false$$otherwise. 316 317 end 318 ----------------------------------------------------------------------------- 319 */ 320 # include <cppad/local/std_set.hpp> 323 324 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 325 /*! 326 \file sparse_hessian.hpp 327 Sparse Hessian driver routine and helper functions. 328 */ 329 // =========================================================================== 330 /*! 331 class used by SparseHessian to hold information 332 so it does not need to be recomputed. 333 */ 335 public: 336 /// Coloring method: "cppad", or "colpack" 337 /// (this field is set by user) 338 std::string color_method; 339 /// row and column indicies for return values 340 /// (some may be reflected by star coloring algorithm) 343 /// indices that sort the user row and col arrays by color 345 /// results of the coloring algorithm 347 348 /// constructor 349 sparse_hessian_work(void) : color_method("cppad.symmetric") 350 { } 351 /// inform CppAD that this information needs to be recomputed 352 void clear(void) 353 { color_method = "cppad.symmetric"; 354 row.clear(); 355 col.clear(); 356 order.clear(); 357 color.clear(); 358 } 359 }; 360 // =========================================================================== 361 /*! 362 Private helper function that does computation for all Sparse Hessian cases. 363 364 \tparam Base 365 is the base type for the recording that is stored in this ADFun<Base object. 366 367 \tparam VectorBase 368 is a simple vector class with elements of type \a Base. 369 370 \tparam VectorSet 371 is a simple vector class with elements of type 372 \c bool or \c std::set<size_t>. 373 374 \tparam VectorSize 375 is sparse_pack or sparse_list. 376 377 \param x [in] 378 is a vector specifing the point at which to compute the Hessian. 379 380 \param w [in] 381 is the weighting vector that defines a scalar valued function by 382 a weighted sum of the components of the vector valued function 383 latex F(x)$$. 384 385 \param sparsity [in] 386 is the sparsity pattern for the Hessian that we are calculating. 387 388 \param user_row [in] 389 is the vector of row indices for the returned Hessian values. 390 391 \param user_col [in] 392 is the vector of columns indices for the returned Hessian values. 393 It must have the same size as user_row. 394 395 \param hes [out] 396 is the vector of Hessian values. 397 It must have the same size as user_row. 398 The return value <code>hes[k]</code> is the second partial of 399 \f$ w^{\rm T} F(x)\f$with respect to the 400 <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. 401 402 \param work 403 This structure contains information that is computed by \c SparseHessianCompute. 404 If the sparsity pattern, \c row vector, or \c col vectors 405 are not the same between calls to \c SparseHessianCompute, 406 \c work.clear() must be called to reinitialize \c work. 407 408 \return 409 Is the number of first order forward sweeps used to compute the 410 requested Hessian values. 411 (This is also equal to the number of second order reverse sweeps.) 412 The total work, not counting the zero order 413 forward sweep, or the time to combine computations, is proportional to this 414 return value. 415 */ 416 template<class Base> 417 template <class VectorBase, class VectorSet, class VectorSize> 419 const VectorBase& x , 420 const VectorBase& w , 421 VectorSet& sparsity , 422 const VectorSize& user_row , 423 const VectorSize& user_col , 424 VectorBase& hes , 425 sparse_hessian_work& work ) 426 { 427 using CppAD::vectorBool; 428 size_t i, k, ell; 429 430 CppAD::vector<size_t>& row(work.row); 431 CppAD::vector<size_t>& col(work.col); 432 CppAD::vector<size_t>& color(work.color); 433 CppAD::vector<size_t>& order(work.order); 434 435 size_t n = Domain(); 436 437 // some values 438 const Base zero(0); 439 const Base one(1); 440 441 // check VectorBase is Simple Vector class with Base type elements 442 CheckSimpleVector<Base, VectorBase>(); 443 444 // number of components of Hessian that are required 445 size_t K = hes.size(); 446 CPPAD_ASSERT_UNKNOWN( size_t( user_row.size() ) == K ); 447 CPPAD_ASSERT_UNKNOWN( size_t( user_col.size() ) == K ); 448 449 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 450 CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); 451 CPPAD_ASSERT_UNKNOWN( row.size() == 0 || row.size() == K ); 452 CPPAD_ASSERT_UNKNOWN( col.size() == 0 || col.size() == K ); 453 454 455 // Point at which we are evaluating the Hessian 456 Forward(0, x); 457 458 // check for case where nothing (except Forward above) to do 459 if( K == 0 ) 460 return 0; 461 462 // Rows of the Hessian (i below) correspond to the forward mode index 463 // and columns (j below) correspond to the reverse mode index. 464 if( color.size() == 0 ) 465 { 466 CPPAD_ASSERT_UNKNOWN( sparsity.n_set() == n ); 467 CPPAD_ASSERT_UNKNOWN( sparsity.end() == n ); 468 469 // copy user rwo and col to work space 470 row.resize(K); 471 col.resize(K); 472 for(k = 0; k < K; k++) 473 { row[k] = user_row[k]; 474 col[k] = user_col[k]; 475 } 476 477 // execute coloring algorithm 478 color.resize(n); 479 if( work.color_method == "cppad.general" ) 480 local::color_general_cppad(sparsity, row, col, color); 481 else if( work.color_method == "cppad.symmetric" ) 482 local::color_symmetric_cppad(sparsity, row, col, color); 483 else if( work.color_method == "colpack.general" ) 484 { 485 # if CPPAD_HAS_COLPACK 486 local::color_general_colpack(sparsity, row, col, color); 487 # else 489 false, 490 "SparseHessian: work.color_method = colpack.general " 491 "and colpack_prefix missing from cmake command line." 492 ); 493 # endif 494 } 495 else if( 496 work.color_method == "colpack.symmetric" || 497 work.color_method == "colpack.star" 498 ) 499 { 500 # if CPPAD_HAS_COLPACK 501 local::color_symmetric_colpack(sparsity, row, col, color); 502 # else 504 false, 505 "SparseHessian: work.color_method is " 506 "colpack.symmetric or colpack.star\n" 507 "and colpack_prefix missing from cmake command line." 508 ); 509 # endif 510 } 511 else 513 false, 514 "SparseHessian: work.color_method is not valid." 515 ); 516 } 517 518 // put sorting indices in color order 519 VectorSize key(K); 520 order.resize(K); 521 for(k = 0; k < K; k++) 522 key[k] = color[ row[k] ]; 523 index_sort(key, order); 524 525 } 526 size_t n_color = 1; 527 for(ell = 0; ell < n; ell++) if( color[ell] < n ) 528 n_color = std::max(n_color, color[ell] + 1); 529 530 // direction vector for calls to forward (rows of the Hessian) 531 VectorBase u(n); 532 533 // location for return values from reverse (columns of the Hessian) 534 VectorBase ddw(2 * n); 535 536 // initialize the return value 537 for(k = 0; k < K; k++) 538 hes[k] = zero; 539 540 // loop over colors 541 # ifndef NDEBUG 542 const std::string& coloring = work.color_method; 543 # endif 544 k = 0; 545 for(ell = 0; ell < n_color; ell++) 546 if( k == K ) 547 { // kludge because colpack returns colors that are not used 548 // (it does not know about the subset corresponding to row, col) 550 coloring == "colpack.general" || 551 coloring == "colpack.symmetic" || 552 coloring == "colpack.star" 553 ); 554 } 555 else if( color[ row[ order[k] ] ] != ell ) 556 { // kludge because colpack returns colors that are not used 557 // (it does not know about the subset corresponding to row, col) 559 coloring == "colpack.general" || 560 coloring == "colpack.symmetic" || 561 coloring == "colpack.star" 562 ); 563 } 564 else 565 { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell ); 566 567 // combine all rows with this color 568 for(i = 0; i < n; i++) 569 { u[i] = zero; 570 if( color[i] == ell ) 571 u[i] = one; 572 } 573 // call forward mode for all these rows at once 574 Forward(1, u); 575 576 // evaluate derivative of w^T * F'(x) * u 577 ddw = Reverse(2, w); 578 579 // set the corresponding components of the result 580 while( k < K && color[ row[ order[k] ] ] == ell ) 581 { hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ]; 582 k++; 583 } 584 } 585 return n_color; 586 } 587 // =========================================================================== 588 // Public Member Functions 589 // =========================================================================== 590 /*! 591 Compute user specified subset of a sparse Hessian. 592 593 The C++ source code corresponding to this operation is 594 \verbatim 595 SparceHessian(x, w, p, row, col, hes, work) 596 \endverbatim 597 598 \tparam Base 599 is the base type for the recording that is stored in this ADFun<Base object. 600 601 \tparam VectorBase 602 is a simple vector class with elements of type \a Base. 603 604 \tparam VectorSet 605 is a simple vector class with elements of type 606 \c bool or \c std::set<size_t>. 607 608 \tparam VectorSize 609 is a simple vector class with elements of type \c size_t. 610 611 \param x [in] 612 is a vector specifing the point at which to compute the Hessian. 613 614 \param w [in] 615 is the weighting vector that defines a scalar valued function by 616 a weighted sum of the components of the vector valued function 617$latex F(x).
618
619 \param p [in]
620 is the sparsity pattern for the Hessian that we are calculating.
621
622 \param row [in]
623 is the vector of row indices for the returned Hessian values.
624
625 \param col [in]
626 is the vector of columns indices for the returned Hessian values.
627 It must have the same size are r.
628
629 \param hes [out]
630 is the vector of Hessian values.
631 It must have the same size are r.
632 The return value <code>hes[k]</code> is the second partial of
633 \f$w^{\rm T} F(x)\f$ with respect to the
634 <code>row[k]</code> and <code>col[k]</code> component of \f$x\f$.
635
636 \param work
637 This structure contains information that is computed by \c SparseHessianCompute.
638 If the sparsity pattern, \c row vector, or \c col vectors
639 are not the same between calls to \c SparseHessian,
640 \c work.clear() must be called to reinitialize \c work.
641
642 \return
643 Is the number of first order forward sweeps used to compute the
644 requested Hessian values.
645 (This is also equal to the number of second order reverse sweeps.)
646 The total work, not counting the zero order
647 forward sweep, or the time to combine computations, is proportional to this
648 return value.
649 */
650 template<class Base>
651 template <class VectorBase, class VectorSet, class VectorSize>
653  const VectorBase& x ,
654  const VectorBase& w ,
655  const VectorSet& p ,
656  const VectorSize& row ,
657  const VectorSize& col ,
658  VectorBase& hes ,
659  sparse_hessian_work& work )
660 {
661  size_t n = Domain();
662  size_t K = hes.size();
663 # ifndef NDEBUG
664  size_t k;
666  size_t(x.size()) == n ,
667  "SparseHessian: size of x not equal domain dimension for f."
668  );
670  size_t(row.size()) == K && size_t(col.size()) == K ,
671  "SparseHessian: either r or c does not have the same size as ehs."
672  );
674  work.color.size() == 0 || work.color.size() == n,
675  "SparseHessian: invalid value in work."
676  );
677  for(k = 0; k < K; k++)
679  row[k] < n,
680  "SparseHessian: invalid value in r."
681  );
683  col[k] < n,
684  "SparseHessian: invalid value in c."
685  );
686  }
687  if( work.color.size() != 0 )
688  for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
689  work.color[j] <= n,
690  "SparseHessian: invalid value in work."
691  );
692 # endif
693  // check for case where there is nothing to compute
694  size_t n_sweep = 0;
695  if( K == 0 )
696  return n_sweep;
697
698  typedef typename VectorSet::value_type Set_type;
699  typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
700  Pattern_type s;
701  if( work.color.size() == 0 )
702  { bool transpose = false;
703  const char* error_msg = "SparseHessian: sparsity pattern"
704  " does not have proper row or column dimension";
705  sparsity_user2internal(s, p, n, n, transpose, error_msg);
706  }
707  n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work);
708  return n_sweep;
709 }
710 /*!
711 Compute a sparse Hessian.
712
713 The C++ source code coresponding to this operation is
714 \verbatim
715  hes = SparseHessian(x, w, p)
716 \endverbatim
717
718
719 \tparam Base
720 is the base type for the recording that is stored in this
722
723 \tparam VectorBase
724 is a simple vector class with elements of the \a Base.
725
726 \tparam VectorSet
727 is a simple vector class with elements of type
728 \c bool or \c std::set<size_t>.
729
730 \param x [in]
731 is a vector specifing the point at which to compute the Hessian.
732
733 \param w [in]
734 The Hessian is computed for a weighted sum of the components
735 of the function corresponding to this ADFun<Base> object.
736 The argument \a w specifies the weights for each component.
737 It must have size equal to the range dimension for this ADFun<Base> object.
738
739 \param p [in]
740 is a sparsity pattern for the Hessian.
741
742 \return
743 Will be a vector of size \c n * n containing the Hessian of
744 at the point specified by \a x
745 (where \c n is the domain dimension for this ADFun<Base> object).
746 */
747 template <class Base>
748 template <class VectorBase, class VectorSet>
750  const VectorBase& x, const VectorBase& w, const VectorSet& p
751 )
752 { size_t i, j, k;
753
754  size_t n = Domain();
755  VectorBase hes(n * n);
756
758  size_t(x.size()) == n,
759  "SparseHessian: size of x not equal domain size for f."
760  );
761
762  typedef typename VectorSet::value_type Set_type;
763  typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
764
765  // initialize the return value as zero
766  Base zero(0);
767  for(i = 0; i < n; i++)
768  for(j = 0; j < n; j++)
769  hes[i * n + j] = zero;
770
771  // arguments to SparseHessianCompute
772  Pattern_type s;
775  sparse_hessian_work work;
776  bool transpose = false;
777  const char* error_msg = "SparseHessian: sparsity pattern"
778  " does not have proper row or column dimension";
779  sparsity_user2internal(s, p, n, n, transpose, error_msg);
780  k = 0;
781  for(i = 0; i < n; i++)
782  { typename Pattern_type::const_iterator itr(s, i);
783  j = *itr;
784  while( j != s.end() )
785  { row.push_back(i);
786  col.push_back(j);
787  k++;
788  j = *(++itr);
789  }
790  }
791  size_t K = k;
792  VectorBase H(K);
793
794  // now we have folded this into the following case
795  SparseHessianCompute(x, w, s, row, col, H, work);
796
797  // now set the non-zero return values
798  for(k = 0; k < K; k++)
799  hes[ row[k] * n + col[k] ] = H[k];
800
801  return hes;
802 }
803 /*!
804 Compute a sparse Hessian
805
806 The C++ source code coresponding to this operation is
807 \verbatim
808  hes = SparseHessian(x, w)
809 \endverbatim
810
811
812 \tparam Base
813 is the base type for the recording that is stored in this
815
816 \tparam VectorBase
817 is a simple vector class with elements of the \a Base.
818
819 \param x [in]
820 is a vector specifing the point at which to compute the Hessian.
821
822 \param w [in]
823 The Hessian is computed for a weighted sum of the components
824 of the function corresponding to this ADFun<Base> object.
825 The argument \a w specifies the weights for each component.
826 It must have size equal to the range dimension for this ADFun<Base> object.
827
828 \return
829 Will be a vector of size \c n * n containing the Hessian of
830 at the point specified by \a x
831 (where \c n is the domain dimension for this ADFun<Base> object).
832 */
833 template <class Base>
834 template <class VectorBase>
835 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
836 { size_t i, j, k;
838
839  size_t m = Range();
840  size_t n = Domain();
841
842  // determine the sparsity pattern p for Hessian of w^T F
843  VectorBool r(n * n);
844  for(j = 0; j < n; j++)
845  { for(k = 0; k < n; k++)
846  r[j * n + k] = false;
847  r[j * n + j] = true;
848  }
849  ForSparseJac(n, r);
850  //
851  VectorBool s(m);
852  for(i = 0; i < m; i++)
853  s[i] = w[i] != 0;
854  VectorBool p = RevSparseHes(n, s);
855
856  // compute sparse Hessian
857  return SparseHessian(x, w, p);
858 }
859
861 # endif
class used by SparseHessian to hold information so it does not need to be recomputed.
Check that exp is true, if not print msg and terminate execution.
void color_general_cppad(const VectorSet &pattern, const VectorSize &row, const VectorSize &col, CppAD::vector< size_t > &color)
Determine which rows of a general sparse matrix can be computed together; i.e., do not have non-zero ...
void clear(void)
inform CppAD that this information needs to be recomputed
row and column indicies for return values (some may be reflected by star coloring algorithm) ...
Coloring algorithm for a symmetric sparse matrix.
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
Coloring algorithm for a general sparse matrix.
void resize(size_t n)
change the number of elements in this vector.
Definition: vector.hpp:399
std::string color_method
Coloring method: &quot;cppad&quot;, or &quot;colpack&quot; (this field is set by user)
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
Two constant standard sets (currently used for concept checking).
void index_sort(const VectorKey &keys, VectorSize &ind)
Compute the indices that sort a vector of keys.
Definition: index_sort.hpp:140
void push_back(const Type &s)
add an element to the back of this vector
Definition: vector.hpp:491
Check that exp is true, if not terminate execution.
indices that sort the user row and col arrays by color
CppAD algorithm for determining which rows of a symmetric sparse matrix can be computed together...
Colpack algorithm for determining which rows of a symmetric sparse matrix can be computed together...
VectorBase SparseHessian(const VectorBase &x, const VectorBase &w)
calculate sparse Hessians
size_t SparseHessianCompute(const VectorBase &x, const VectorBase &w, VectorSet &sparsity, const VectorSize &row, const VectorSize &col, VectorBase &hes, sparse_hessian_work &work)
Private helper function that does computation for all Sparse Hessian cases.