CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
rev_sparse_hes.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_REV_SPARSE_HES_HPP
2 # define CPPAD_CORE_REV_SPARSE_HES_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
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>
193 # include <cppad/local/pod_vector.hpp>
194 # include <cppad/local/std_set.hpp>
195 
196 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
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++)
269  { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
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++)
299  { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
300 
301  // ind_taddr_[j] is operator taddr for j-th independent variable
302  CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
303  CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );
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  );
397  CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
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++)
421  { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
422  CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
423  CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );
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
549 in this ADFUN object.
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++)
572  { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ )
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++)
598  { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
599 
600  // ind_taddr_[j] is operator taddr for j-th independent variable
601  CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
602  CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );
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 
619 } // END_CPPAD_NAMESPACE
620 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
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
#define CPPAD_ASSERT_UNKNOWN(exp)
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.