CppAD: A C++ Algorithmic Differentiation Package  20171217
rev_sparse_hes.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 RevSparseHes$$17 spell 18 std 19 VecAD 20 Jacobian 21 Jac 22 Hessian 23 Hes 24 const 25 Bool 26 Dep 27 proportional 28 var 29 cpp 30$$ 31 32$section Hessian Sparsity Pattern: Reverse Mode$$33 34 head Syntax$$
35 $icode%h% = %f%.RevSparseHes(%q%, %s%) 36 %$$37 icode%h% = %f%.RevSparseHes(%q%, %s%, %transpose%)%$$ 38 39$head Purpose$$40 We use latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
41 $cref/AD function/glossary/AD Function/$$corresponding to icode f$$. 42 For a fixed matrix$latex R \in \B{R}^{n \times q}$$43 and a fixed vector latex S \in \B{R}^{1 \times m}$$,
44 we define
45 $latex $46 \begin{array}{rcl} 47 H(x) 48 & = & \partial_x \left[ \partial_u S * F[ x + R * u ] \right]_{u=0} 49 \\ 50 & = & R^\R{T} * (S * F)^{(2)} ( x ) 51 \\ 52 H(x)^\R{T} 53 & = & (S * F)^{(2)} ( x ) * R 54 \end{array} 55$ $$56 Given a 57 cref/sparsity pattern/glossary/Sparsity Pattern/$$ 58 for the matrix$latex R$$and the vector latex S$$,
59 $code RevSparseHes$$returns a sparsity pattern for the latex H(x)$$. 60 61$head f$$62 The object icode f$$ has prototype
63 $codei% 64 const ADFun<%Base%> %f% 65 %$$66 67 head x$$ 68 If the operation sequence in$icode f$$is 69 cref/independent/glossary/Operation/Independent/$$ of
70 the independent variables in $latex x \in B^n$$, 71 the sparsity pattern is valid for all values of 72 (even if it has cref CondExp$$ or$cref VecAD$$operations). 73 74 head q$$
75 The argument $icode q$$has prototype 76 codei% 77 size_t %q% 78 %$$ 79 It specifies the number of columns in$latex R \in \B{R}^{n \times q}$$80 and the number of rows in latex H(x) \in \B{R}^{q \times n}$$.
81 It must be the same value as in the previous $cref ForSparseJac$$call 82 codei% 83 %f%.ForSparseJac(%q%, %r%, %r_transpose%) 84 %$$ 85 Note that if$icode r_transpose$$is true, icode r$$ in the call above
86 corresponding to $latex R^\R{T} \in \B{R}^{q \times n}$$87 88 head transpose$$ 89 The argument$icode transpose$$has prototype 90 codei% 91 bool %transpose% 92 %$$
93 The default value $code false$$is used when icode transpose$$ is not present. 94 95 96$head r$$97 The matrix latex R$$ is specified by the previous call
98 $codei% 99 %f%.ForSparseJac(%q%, %r%, %transpose%) 100 %$$101 see cref/r/ForSparseJac/r/$$. 102 The type of the elements of 103$cref/VectorSet/RevSparseHes/VectorSet/$$must be the 104 same as the type of the elements of icode r$$.
105
106 $head s$$107 The argument icode s$$ has prototype 108$codei%
109  const %VectorSet%& %s%
110 %$$111 (see cref/VectorSet/RevSparseHes/VectorSet/$$ below)
112 If it has elements of type $code bool$$, 113 its size is latex m$$. 114 If it has elements of type$code std::set<size_t>$$, 115 its size is one and all the elements of icode%s%[0]%$$
116 are between zero and $latex m - 1$$. 117 It specifies a 118 cref/sparsity pattern/glossary/Sparsity Pattern/$$ 119 for the vector$icode S$$. 120 121 head h$$
122 The result $icode h$$has prototype 123 codei% 124 %VectorSet%& %h% 125 %$$ 126 (see$cref/VectorSet/RevSparseHes/VectorSet/$$below). 127 128 subhead transpose false$$
129 If $icode h$$has elements of type code bool$$, 130 its size is$latex q * n$$. 131 If it has elements of type code std::set<size_t>$$,
132 its size is $latex q$$and all the set elements are between 133 zero and icode%n%-1%$$ inclusive. 134 It specifies a 135$cref/sparsity pattern/glossary/Sparsity Pattern/$$136 for the matrix latex H(x)$$.
137
138 $subhead transpose true$$139 If icode h$$ has elements of type$code bool$$, 140 its size is latex n * q$$.
141 If it has elements of type $code std::set<size_t>$$, 142 its size is latex n$$ and all the set elements are between 143 zero and$icode%q%-1%$$inclusive. 144 It specifies a 145 cref/sparsity pattern/glossary/Sparsity Pattern/$$
146 for the matrix $latex H(x)^\R{T}$$. 147 148 head VectorSet$$ 149 The type$icode VectorSet$$must be a cref SimpleVector$$ class with
150 $cref/elements of type/SimpleVector/Elements of Specified Type/$$151 code bool$$ or$code std::set<size_t>$$; 152 see cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
153 of the difference.
154 The type of the elements of
155 $cref/VectorSet/RevSparseHes/VectorSet/$$must be the 156 same as the type of the elements of icode r$$. 157 158$head Entire Sparsity Pattern$$159 Suppose that latex q = n$$ and
160 $latex R \in \B{R}^{n \times n}$$is the latex n \times n$$ identity matrix. 161 Further suppose that the$latex S$$is the th k$$
162 $cref/elementary vector/glossary/Elementary Vector/$$; i.e. 163 latex $164 S_j = \left\{ \begin{array}{ll} 165 1 & {\rm if} \; j = k 166 \\ 167 0 & {\rm otherwise} 168 \end{array} \right. 169$$$ 170 In this case, 171 the corresponding value$icode h$$is a 172 sparsity pattern for the Hessian matrix 173 latex F_k^{(2)} (x) \in \B{R}^{n \times n}$$.
174
175 $head Example$$176 children% 177 example/sparse/rev_sparse_hes.cpp 178 %example/sparse/sparsity_sub.cpp 179 %$$ 180 The file 181$cref rev_sparse_hes.cpp$$182 contains an example and test of this operation. 183 It returns true if it succeeds and false otherwise. 184 The file 185 cref/sparsity_sub.cpp/sparsity_sub.cpp/RevSparseHes/$$
186 contains an example and test of using $code RevSparseHes$$187 to compute the sparsity pattern for a subset of the Hessian. 188 189$end
190 -----------------------------------------------------------------------------
191 */
192 # include <algorithm>
195
197 /*!
198 \file rev_sparse_hes.hpp
199 Reverse mode Hessian sparsity patterns.
200 */
201 // ===========================================================================
202 // RevSparseHesCase
203 /*!
204 Private helper function for RevSparseHes(q, s) bool sparsity.
205
206 All of the description in the public member function RevSparseHes(q, s)
207 applies.
208
209 \param set_type
210 is a \c bool value. This argument is used to dispatch to the proper source
211 code depending on the vlaue of \c VectorSet::value_type.
212
213 \param transpose
214 See \c RevSparseHes(q, s).
215
216 \param q
217 See \c RevSparseHes(q, s).
218
219 \param s
220 See \c RevSparseHes(q, s).
221
222 \param h
223 is the return value for the corresponging call to \c RevSparseJac(q, s).
224 */
225 template <class Base>
226 template <class VectorSet>
228  bool set_type ,
229  bool transpose ,
230  size_t q ,
231  const VectorSet& s ,
232  VectorSet& h )
233 { size_t n = Domain();
234  size_t m = Range();
235  //
236  h.resize(q * n );
237
239  for_jac_sparse_pack_.n_set() > 0,
240  "RevSparseHes: previous stored call to ForSparseJac did not "
241  "use bool for the elements of r."
242  );
243  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );
244  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == num_var_tape_ );
245  //
246  // temporary indices
247  size_t i, j;
248
249  // check Vector is Simple VectorSet class with bool elements
250  CheckSimpleVector<bool, VectorSet>();
251
253  q == for_jac_sparse_pack_.end(),
254  "RevSparseHes: q is not equal to its value\n"
255  "in the previous call to ForSparseJac with this ADFun object."
256  );
258  size_t(s.size()) == m,
259  "RevSparseHes: size of s is not equal to\n"
260  "range dimension for ADFun object."
261  );
262
263  // Array that will hold reverse Jacobian dependency flag.
264  // Initialize as true for the dependent variables.
265  local::pod_vector<bool> RevJac(num_var_tape_);
266  for(i = 0; i < num_var_tape_; i++)
267  RevJac[i] = false;
268  for(i = 0; i < m; i++)
270  RevJac[ dep_taddr_[i] ] = s[i];
271  }
272
273  // vector of sets that will hold reverse Hessain values
274  local::sparse_pack rev_hes_pattern;
275  rev_hes_pattern.resize(num_var_tape_, q);
276
277  // compute the Hessian sparsity patterns
279  &play_,
280  n,
281  num_var_tape_,
282  for_jac_sparse_pack_,
283  RevJac.data(),
284  rev_hes_pattern
285  );
286
287  // return values corresponding to independent variables
288  CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == n * q );
289  for(j = 0; j < n; j++)
290  { for(i = 0; i < q; i++)
291  { if( transpose )
292  h[ j * q + i ] = false;
293  else h[ i * n + j ] = false;
294  }
295  }
296
297  // j is index corresponding to reverse mode partial
298  for(j = 0; j < n; j++)
300
304
305  // extract the result from rev_hes_pattern
306  CPPAD_ASSERT_UNKNOWN( rev_hes_pattern.end() == q );
307  local::sparse_pack::const_iterator itr(rev_hes_pattern, j + 1);
308  i = *itr;
309  while( i < q )
310  { if( transpose )
311  h[ j * q + i ] = true;
312  else h[ i * n + j ] = true;
313  i = *(++itr);
314  }
315  }
316 }
317 /*!
318 Private helper function for RevSparseHes(q, s) set sparsity.
319
320 All of the description in the public member function RevSparseHes(q, s)
321 applies.
322
323 \param set_type
324 is a \c std::set<size_t> value.
325 This argument is used to dispatch to the proper source
326 code depending on the vlaue of \c VectorSet::value_type.
327
328 \param transpose
329 See \c RevSparseHes(q, s).
330
331 \param q
332 See \c RevSparseHes(q, s).
333
334 \param s
335 See \c RevSparseHes(q, s).
336
337 \param h
338 is the return value for the corresponging call to \c RevSparseJac(q, s).
339 */
340 template <class Base>
341 template <class VectorSet>
343  const std::set<size_t>& set_type ,
344  bool transpose ,
345  size_t q ,
346  const VectorSet& s ,
347  VectorSet& h )
348 { size_t n = Domain();
349 # ifndef NDEBUG
350  size_t m = Range();
351 # endif
352  //
353  if( transpose )
354  h.resize(n);
355  else h.resize(q);
356
358  for_jac_sparse_set_.n_set() > 0,
359  "RevSparseHes: previous stored call to ForSparseJac did not "
360  "use std::set<size_t> for the elements of r."
361  );
362  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
363  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == num_var_tape_ );
364  //
365  // temporary indices
366  size_t i, j;
367  std::set<size_t>::const_iterator itr_1;
368
369  // check VectorSet is Simple Vector class with sets for elements
370  CheckSimpleVector<std::set<size_t>, VectorSet>(
371  local::one_element_std_set<size_t>(), local::two_element_std_set<size_t>()
372  );
373
375  q == for_jac_sparse_set_.end(),
376  "RevSparseHes: q is not equal to its value\n"
377  "in the previous call to ForSparseJac with this ADFun object."
378  );
380  s.size() == 1,
381  "RevSparseHes: size of s is not equal to one."
382  );
383
384  // Array that will hold reverse Jacobian dependency flag.
385  // Initialize as true for the dependent variables.
386  local::pod_vector<bool> RevJac(num_var_tape_);
387  for(i = 0; i < num_var_tape_; i++)
388  RevJac[i] = false;
389  itr_1 = s[0].begin();
390  while( itr_1 != s[0].end() )
391  { i = *itr_1++;
393  i < m,
394  "RevSparseHes: an element of the set s[0] has value "
395  "greater than or equal m"
396  );
398  RevJac[ dep_taddr_[i] ] = true;
399  }
400
401
402  // vector of sets that will hold reverse Hessain values
403  local::sparse_list rev_hes_pattern;
404  rev_hes_pattern.resize(num_var_tape_, q);
405
406  // compute the Hessian sparsity patterns
408  &play_,
409  n,
410  num_var_tape_,
411  for_jac_sparse_set_,
412  RevJac.data(),
413  rev_hes_pattern
414  );
415
416  // return values corresponding to independent variables
417  // j is index corresponding to reverse mode partial
418  CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == q || transpose );
419  CPPAD_ASSERT_UNKNOWN( size_t(h.size()) == n || ! transpose );
420  for(j = 0; j < n; j++)
424
425  // extract the result from rev_hes_pattern
426  // and add corresponding elements to result sets in h
427  CPPAD_ASSERT_UNKNOWN( rev_hes_pattern.end() == q );
428  local::sparse_list::const_iterator itr_2(rev_hes_pattern, j+1);
429  i = *itr_2;
430  while( i < q )
431  { if( transpose )
432  h[j].insert(i);
433  else h[i].insert(j);
434  i = *(++itr_2);
435  }
436  }
437 }
438
439 // ===========================================================================
440 // RevSparseHes
441
442 /*!
443 User API for Hessian sparsity patterns using reverse mode.
444
445 The C++ source code corresponding to this operation is
446 \verbatim
447  h = f.RevSparseHes(q, r)
448 \endverbatim
449
450 \tparam Base
451 is the base type for this recording.
452
453 \tparam VectorSet
454 is a simple vector with elements of type \c bool
455 or \c std::set<size_t>.
456
457 \param transpose
458 is true (false) if \c is is equal to \f$H(x) \f$ (\f$H(x)^T \f$)
459 where
460 \f[
461  H(x) = R^T (S * F)^{(2)} (x)
462 \f]
463 where \f$F \f$ is the function corresponding to the operation sequence
464 and \a x is any argument value.
465
466 \param q
467 is the value of \a q in the
468 by the previous call of the form
469 \verbatim
470  f.ForSparseJac(q, r, packed)
471 \endverbatim
472 The value \c r in this call is a sparsity pattern for the matrix \f$R \f$.
473 The type of the element of \c r for the previous call to \c ForSparseJac
474 must be the same as the type of the elements of \c s.
475
476 \param s
477 is a vector with size \c m that specifies the sparsity pattern
478 for the vector \f$S \f$,
479 where \c m is the number of dependent variables
480 corresponding to the operation sequence stored in \a play.
481
482 \return
483 If \c transpose is false (true),
484 the return vector is a sparsity pattern for \f$H(x) \f$ (\f$H(x)^T \f$).
485 \f[
486  H(x) = R^T ( S * F)^{(2)} (x)
487 \f]
488 where \f$F \f$ is the function corresponding to the operation sequence
489 and \a x is any argument value.
490 */
491
492 template <class Base>
493 template <class VectorSet>
495  size_t q, const VectorSet& s, bool transpose
496 )
497 { VectorSet h;
498  typedef typename VectorSet::value_type Set_type;
499
500  // Should check to make sure q is same as in previous call to
501  // forward sparse Jacobian.
502  RevSparseHesCase(
503  Set_type() ,
504  transpose ,
505  q ,
506  s ,
507  h
508  );
509
510  return h;
511 }
512 // ===========================================================================
513 // RevSparseHesCheckpoint
514 /*!
515 Hessian sparsity patterns calculation used by checkpoint functions.
516
517 \tparam Base
518 is the base type for this recording.
519
520 \param transpose
521 is true (false) h is equal to \f$H(x) \f$ (\f$H(x)^T \f$)
522 where
523 \f[
524  H(x) = R^T (S * F)^{(2)} (x)
525 \f]
526 where \f$F \f$ is the function corresponding to the operation sequence
527 and \f$x \f$ is any argument value.
528
529 \param q
530 is the value of q in the by the previous call of the form
531 \verbatim
532  f.ForSparseJac(q, r)
533 \endverbatim
534 The value r in this call is a sparsity pattern for the matrix \f$R \f$.
535
536 \param s
537 is a vector with size m that specifies the sparsity pattern
538 for the vector \f$S \f$,
539 where m is the number of dependent variables
540 corresponding to the operation sequence stored in play_.
541
542 \param h
543 The input size and elements of h do not matter.
544 On output, h is the sparsity pattern for the matrix \f$H(x) \f$
545 or \f$H(x)^T \f$ depending on transpose.
546
547 \par Assumptions
548 The forward jacobian sparsity pattern must be currently stored
550 */
551 template <class Base>
553  size_t q ,
554  vector<bool>& s ,
555  bool transpose ,
556  local::sparse_list& h )
557 { size_t n = Domain();
558  size_t m = Range();
559
560  // checkpoint functions should get this right
561  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
562  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == num_var_tape_ );
563  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.end() == q );
564  CPPAD_ASSERT_UNKNOWN( s.size() == m );
565
566  // Array that holds the reverse Jacobiain dependcy flags.
567  // Initialize as true for dependent variables, flase for others.
568  local::pod_vector<bool> RevJac(num_var_tape_);
569  for(size_t i = 0; i < num_var_tape_; i++)
570  RevJac[i] = false;
571  for(size_t i = 0; i < m; i++)
573  RevJac[ dep_taddr_[i] ] = s[i];
574  }
575
576  // holds reverse Hessian sparsity pattern for all variables
577  local::sparse_list rev_hes_pattern;
578  rev_hes_pattern.resize(num_var_tape_, q);
579
580  // compute Hessian sparsity pattern for all variables
582  &play_,
583  n,
584  num_var_tape_,
585  for_jac_sparse_set_,
586  RevJac.data(),
587  rev_hes_pattern
588  );
589
590  // dimension the return value
591  if( transpose )
592  h.resize(n, q);
593  else
594  h.resize(q, n);
595
596  // j is index corresponding to reverse mode partial
597  for(size_t j = 0; j < n; j++)
599
603
604  // extract the result from rev_hes_pattern
605  CPPAD_ASSERT_UNKNOWN( rev_hes_pattern.end() == q );
606  local::sparse_list::const_iterator itr(rev_hes_pattern, j + 1);
607  size_t i = *itr;
608  while( i < q )
609  { if( transpose )
610  h.post_element(j, i);
611  else h.post_element(i, j);
612  i = *(++itr);
613  }
614  }
615  for(size_t i = 0; i < h.n_set(); ++i)
616  h.process_post(i);
617 }
618
620 # endif
Check that exp is true, if not print msg and terminate execution.
size_t end(void) const
Fetch end for this vector of sets object.
Vector of sets of positive integers, each set stored as a singly linked list.
Definition: sparse_list.hpp:35
void RevSparseHesCheckpoint(size_t q, vector< bool > &s, bool transpose, local::sparse_list &h)
Hessian sparsity patterns calculation used by checkpoint functions.
size_t end(void) const
Fetch end for this vector of sets object.
File used to define pod_vector class.
void rev_hes_sweep(const local::player< Base > *play, size_t n, size_t numvar, const Vector_set &for_jac_sparse, bool *RevJac, Vector_set &rev_hes_sparse)
Given the forward Jacobian sparsity pattern for all the variables, and the reverse Jacobian sparsity ...
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
Vector of sets of postivie integers, each set stored as a packed boolean array.
Definition: sparse_pack.hpp:33
VectorSet RevSparseHes(size_t q, const VectorSet &s, bool transpose=false)
User API for Hessian sparsity patterns using reverse mode.
void post_element(size_t i, size_t element)
Post an element for delayed addition to a set.
Two constant standard sets (currently used for concept checking).
Type * data(void)
current data pointer is no longer valid after any of the following: extend, erase, operator=, and ~pod_vector. Take extreem care when using this function.
Definition: pod_vector.hpp:89
Check that exp is true, if not terminate execution.
void resize(size_t n_set, size_t end)
Start a new vector of sets.
cons_iterator for one set of positive integers in a sparse_pack object.
void process_post(size_t i)
process post entries for a specific set.
size_t n_set(void) const
Fetch n_set for vector of sets object.
cons_iterator for one set of positive integers in a sparse_list object.
void RevSparseHesCase(bool set_type, bool transpose, size_t q, const VectorSet &s, VectorSet &h)
Private helper function for RevSparseHes(q, s) bool sparsity.
Scalar value_type
void resize(size_t n_set, size_t end)
Change number of sets, set end, and initialize all sets as empty.