CppAD: A C++ Algorithmic Differentiation Package  20171217
for_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 ForSparseHes$$17 spell 18 Andrea Walther 19 std 20 VecAD 21 Jacobian 22 Jac 23 Hessian 24 Hes 25 const 26 Bool 27 Dep 28 proportional 29 var 30 cpp 31$$ 32 33$section Hessian Sparsity Pattern: Forward Mode$$34 35 head Syntax$$
36 $icode%h% = %f%.ForSparseHes(%r%, %s%) 37 %$$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 we define 43 latex $44 \begin{array}{rcl} 45 H(x) 46 & = & \partial_x \left[ \partial_u S \cdot F[ x + R \cdot u ] \right]_{u=0} 47 \\ 48 & = & R^\R{T} \cdot (S \cdot F)^{(2)} ( x ) \cdot R 49 \end{array} 50$$$ 51 Where$latex R \in \B{R}^{n \times n}$$is a diagonal matrix 52 and latex S \in \B{R}^{1 \times m}$$ is a row vector.
53 Given a
54 $cref/sparsity pattern/glossary/Sparsity Pattern/$$55 for the diagonal of latex R$$ and the vector$latex S$$, 56 code ForSparseHes$$ returns a sparsity pattern for the $latex H(x)$$. 57 58 head f$$ 59 The object$icode f$$has prototype 60 codei% 61 const ADFun<%Base%> %f% 62 %$$
63
64 $head x$$65 If the operation sequence in icode f$$ is 66$cref/independent/glossary/Operation/Independent/$$of 67 the independent variables in latex x \in B^n$$,
68 the sparsity pattern is valid for all values of
69 (even if it has $cref CondExp$$or cref VecAD$$ operations). 70 71$head r$$72 The argument icode r$$ has prototype
73 $codei% 74 const %VectorSet%& %r% 75 %$$76 (see cref/VectorSet/ForSparseHes/VectorSet/$$ below) 77 If it has elements of type$code bool$$, 78 its size is latex n$$.
79 If it has elements of type $code std::set<size_t>$$, 80 its size is one and all the elements of icode%s%[0]%$$ 81 are between zero and$latex n - 1$$. 82 It specifies a 83 cref/sparsity pattern/glossary/Sparsity Pattern/$$
84 for the diagonal of $latex R$$. 85 The fewer non-zero elements in this sparsity pattern, 86 the faster the calculation should be and the more sparse 87 latex H(x)$$ should be. 88 89$head s$$90 The argument icode s$$ has prototype
91 $codei% 92 const %VectorSet%& %s% 93 %$$94 (see cref/VectorSet/ForSparseHes/VectorSet/$$ below) 95 If it has elements of type$code bool$$, 96 its size is latex m$$.
97 If it has elements of type $code std::set<size_t>$$, 98 its size is one and all the elements of icode%s%[0]%$$ 99 are between zero and$latex m - 1$$. 100 It specifies a 101 cref/sparsity pattern/glossary/Sparsity Pattern/$$
102 for the vector $icode S$$. 103 The fewer non-zero elements in this sparsity pattern, 104 the faster the calculation should be and the more sparse 105 latex H(x)$$ should be. 106 107$head h$$108 The result icode h$$ has prototype
109 $codei% 110 %VectorSet%& %h% 111 %$$112 (see cref/VectorSet/ForSparseHes/VectorSet/$$ below). 113 If$icode h$$has elements of type code bool$$,
114 its size is $latex n * n$$. 115 If it has elements of type code std::set<size_t>$$, 116 its size is$latex n$$and all the set elements are between 117 zero and icode%n%-1%$$ inclusive.
118 It specifies a
119 $cref/sparsity pattern/glossary/Sparsity Pattern/$$120 for the matrix latex H(x)$$. 121 122$head VectorSet$$123 The type icode VectorSet$$ must be a $cref SimpleVector$$class with 124 cref/elements of type/SimpleVector/Elements of Specified Type/$$ 125$code bool$$or code std::set<size_t>$$;
126 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$for a discussion 127 of the difference. 128 The type of the elements of 129 cref/VectorSet/ForSparseHes/VectorSet/$$ must be the 130 same as the type of the elements of$icode r$$. 131 132 head Algorithm$$
133 See Algorithm II in
134 $italic Computing sparse Hessians with automatic differentiation$$135 by Andrea Walther. 136 Note that icode s$$ provides the information so that 137 'dead ends' are not included in the sparsity pattern. 138 139$head Example$$140 children% 141 example/sparse/for_sparse_hes.cpp 142 %$$
143 The file
144 $cref for_sparse_hes.cpp$$145 contains an example and test of this operation. 146 It returns true if it succeeds and false otherwise. 147 148$end
149 -----------------------------------------------------------------------------
150 */
151 # include <algorithm>
154
156 /*!
157 \file for_sparse_hes.hpp
158 Forward mode Hessian sparsity patterns.
159 */
160 // ===========================================================================
161 // ForSparseHesCase
162 /*!
163 Private helper function for ForSparseHes(q, s) bool sparsity.
164
165 All of the description in the public member function ForSparseHes(q, s)
166 applies.
167
168 \param set_type
169 is a \c bool value. This argument is used to dispatch to the proper source
170 code depending on the vlaue of \c VectorSet::value_type.
171
172 \param r
173 See \c ForSparseHes(r, s).
174
175 \param s
176 See \c ForSparseHes(r, s).
177
178 \param h
179 is the return value for the corresponging call to \c ForSparseJac(q, s).
180 */
181 template <class Base>
182 template <class VectorSet>
184  bool set_type ,
185  const VectorSet& r ,
186  const VectorSet& s ,
187  VectorSet& h )
188 { size_t n = Domain();
189  size_t m = Range();
190  //
191  // check Vector is Simple VectorSet class with bool elements
192  CheckSimpleVector<bool, VectorSet>();
193  //
195  size_t(r.size()) == n,
196  "ForSparseHes: size of r is not equal to\n"
197  "domain dimension for ADFun object."
198  );
200  size_t(s.size()) == m,
201  "ForSparseHes: size of s is not equal to\n"
202  "range dimension for ADFun object."
203  );
204  //
205  // sparsity pattern corresponding to r
206  local::sparse_pack for_jac_pattern;
207  for_jac_pattern.resize(num_var_tape_, n + 1);
208  for(size_t i = 0; i < n; i++)
212  //
214  if( r[i] )
216  }
217  // compute forward Jacobiain sparsity pattern
218  bool dependency = false;
220  &play_,
221  dependency,
222  n,
223  num_var_tape_,
224  for_jac_pattern
225  );
226  // sparsity pattern correspnding to s
227  local::sparse_pack rev_jac_pattern;
228  rev_jac_pattern.resize(num_var_tape_, 1);
229  for(size_t i = 0; i < m; i++)
231  //
233  if( s[i] )
235  }
236  // compute reverse sparsity pattern for dependency analysis
237  // (note that we are only want non-zero derivatives not true dependency)
239  &play_,
240  dependency,
241  n,
242  num_var_tape_,
243  rev_jac_pattern
244  );
245  // vector of sets that will hold the forward Hessain values
246  local::sparse_pack for_hes_pattern;
247  for_hes_pattern.resize(n+1, n+1);
248  //
249  // compute the Hessian sparsity patterns
251  &play_,
252  n,
253  num_var_tape_,
254  for_jac_pattern,
255  rev_jac_pattern,
256  for_hes_pattern
257  );
258  // initialize return values corresponding to independent variables
259  h.resize(n * n);
260  for(size_t i = 0; i < n; i++)
261  { for(size_t j = 0; j < n; j++)
262  h[ i * n + j ] = false;
263  }
264  // copy to result pattern
265  CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == n+1 );
266  for(size_t i = 0; i < n; i++)
270
271  // extract the result from for_hes_pattern
273  size_t j = *itr;
274  while( j < for_hes_pattern.end() )
275  { CPPAD_ASSERT_UNKNOWN( 0 < j )
276  h[ i * n + (j-1) ] = true;
277  j = *(++itr);
278  }
279  }
280 }
281 /*!
282 Private helper function for ForSparseHes(q, s) set sparsity.
283
284 All of the description in the public member function ForSparseHes(q, s)
285 applies.
286
287 \param set_type
288 is a \c std::set<size_t> value.
289 This argument is used to dispatch to the proper source
290 code depending on the vlaue of \c VectorSet::value_type.
291
292 \param r
293 See \c ForSparseHes(r, s).
294
295 \param s
296 See \c ForSparseHes(q, s).
297
298 \param h
299 is the return value for the corresponging call to \c ForSparseJac(q, s).
300 */
301 template <class Base>
302 template <class VectorSet>
304  const std::set<size_t>& set_type ,
305  const VectorSet& r ,
306  const VectorSet& s ,
307  VectorSet& h )
308 { size_t n = Domain();
309 # ifndef NDEBUG
310  size_t m = Range();
311 # endif
312  std::set<size_t>::const_iterator itr_1;
313  //
314  // check VectorSet is Simple Vector class with sets for elements
315  CheckSimpleVector<std::set<size_t>, VectorSet>(
316  local::one_element_std_set<size_t>(), local::two_element_std_set<size_t>()
317  );
319  r.size() == 1,
320  "ForSparseHes: size of s is not equal to one."
321  );
323  s.size() == 1,
324  "ForSparseHes: size of s is not equal to one."
325  );
326  //
327  // sparsity pattern corresponding to r
328  local::sparse_list for_jac_pattern;
329  for_jac_pattern.resize(num_var_tape_, n + 1);
330  itr_1 = r[0].begin();
331  while( itr_1 != r[0].end() )
332  { size_t i = *itr_1++;
336  //
339  }
340  // compute forward Jacobiain sparsity pattern
341  bool dependency = false;
343  &play_,
344  dependency,
345  n,
346  num_var_tape_,
347  for_jac_pattern
348  );
349  // sparsity pattern correspnding to s
350  local::sparse_list rev_jac_pattern;
351  rev_jac_pattern.resize(num_var_tape_, 1);
352  itr_1 = s[0].begin();
353  while( itr_1 != s[0].end() )
354  { size_t i = *itr_1++;
356  i < m,
357  "ForSparseHes: an element of the set s[0] has value "
358  "greater than or equal m"
359  );
361  //
364  }
365  //
366  // compute reverse sparsity pattern for dependency analysis
367  // (note that we are only want non-zero derivatives not true dependency)
369  &play_,
370  dependency,
371  n,
372  num_var_tape_,
373  rev_jac_pattern
374  );
375  //
376  // vector of sets that will hold reverse Hessain values
377  local::sparse_list for_hes_pattern;
378  for_hes_pattern.resize(n+1, n+1);
379  //
380  // compute the Hessian sparsity patterns
382  &play_,
383  n,
384  num_var_tape_,
385  for_jac_pattern,
386  rev_jac_pattern,
387  for_hes_pattern
388  );
389  // return values corresponding to independent variables
390  // j is index corresponding to reverse mode partial
391  h.resize(n);
392  CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == n+1 );
393  for(size_t i = 0; i < n; i++)
396
397  // extract the result from for_hes_pattern
399  size_t j = *itr_2;
400  while( j < for_hes_pattern.end() )
401  { CPPAD_ASSERT_UNKNOWN( 0 < j )
402  h[i].insert(j-1);
403  j = *(++itr_2);
404  }
405  }
406 }
407
408 // ===========================================================================
409 // ForSparseHes
410
411 /*!
412 User API for Hessian sparsity patterns using reverse mode.
413
414 The C++ source code corresponding to this operation is
415 \verbatim
416  h = f.ForSparseHes(q, r)
417 \endverbatim
418
419 \tparam Base
420 is the base type for this recording.
421
422 \tparam VectorSet
423 is a simple vector with elements of type \c bool
424 or \c std::set<size_t>.
425
426 \param r
427 is a vector with size \c n that specifies the sparsity pattern
428 for the diagonal of the matrix \f$R \f$,
429 where \c n is the number of independent variables
430 corresponding to the operation sequence stored in \a play.
431
432 \param s
433 is a vector with size \c m that specifies the sparsity pattern
434 for the vector \f$S \f$,
435 where \c m is the number of dependent variables
436 corresponding to the operation sequence stored in \a play.
437
438 \return
439 The return vector is a sparsity pattern for \f$H(x) \f$
440 \f[
441  H(x) = R^T ( S * F)^{(2)} (x) R
442 \f]
443 where \f$F \f$ is the function corresponding to the operation sequence
444 and \a x is any argument value.
445 */
446
447 template <class Base>
448 template <class VectorSet>
450  const VectorSet& r, const VectorSet& s
451 )
452 { VectorSet h;
453  typedef typename VectorSet::value_type Set_type;
454
455  // Should check to make sure q is same as in previous call to
456  // forward sparse Jacobian.
457  ForSparseHesCase(
458  Set_type() ,
459  r ,
460  s ,
461  h
462  );
463
464  return h;
465 }
466 // ===========================================================================
467 // ForSparseHesCheckpoint
468 /*!
469 Hessian sparsity patterns calculation used by checkpoint functions.
470
471 \tparam Base
472 is the base type for this recording.
473
474 \param r
475 is a vector with size n that specifies the sparsity pattern
476 for the diagonal of \f$R \f$,
477 where n is the number of independent variables
478 corresponding to the operation sequence stored in play_.
479
480 \param s
481 is a vector with size m that specifies the sparsity pattern
482 for the vector \f$S \f$,
483 where m is the number of dependent variables
484 corresponding to the operation sequence stored in play_.
485
486 \param h
487 The input size and elements of h do not matter.
488 On output, h is the sparsity pattern for the matrix \f$H(x) R \f$.
489
490 \par Assumptions
491 The forward jacobian sparsity pattern must be currently stored
493 */
494
495 // The checkpoint class is not yet using forward sparse Hessians.
497 template <class Base>
499  vector<bool>& r ,
500  vector<bool>& s ,
501  local::sparse_list& h )
502 {
503  size_t n = Domain();
504  size_t m = Range();
505
506  // checkpoint functions should get this right
507  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
508  CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );
509  CPPAD_ASSERT_UNKNOWN( s.size() == m );
510
511  // Array that holds the reverse Jacobiain dependcy flags.
512  // Initialize as true for dependent variables, flase for others.
513  local::pod_vector<bool> RevJac(num_var_tape_);
514  for(size_t i = 0; i < num_var_tape_; i++)
515  RevJac[i] = false;
516  for(size_t i = 0; i < m; i++)
518  RevJac[ dep_taddr_[i] ] = s[i];
519  }
520
521  // holds forward Hessian sparsity pattern for all variables
522  local::sparse_list for_hes_pattern;
523  for_hes_pattern.resize(n+1, n+1);
524
525  // compute Hessian sparsity pattern for all variables
526  local::for_hes_sweep(
527  &play_,
528  n,
529  num_var_tape_,
530  for_jac_sparse_set_,
531  RevJac.data(),
532  for_hes_pattern
533  );
534
535  // dimension the return value
536  if( transpose )
537  h.resize(n, n);
538  else
539  h.resize(n, n);
540
541  // j is index corresponding to reverse mode partial
542  for(size_t j = 0; j < n; j++)
544
548
549  // extract the result from for_hes_pattern
550  CPPAD_ASSERT_UNKNOWN( for_hes_pattern.end() == q );
551  local::sparse_list::const_iterator itr(for_hes_pattern, .j + 1);
552  size_t i = *itr;
553  while( i < q )
554  { if( transpose )
555  h.post_element(j, i);
556  else h.post_element(i, j);
557  i = *(++itr);
558  }
559  }
560  // process posts
561  for(size_t i = 0; i < n; ++i)
562  h.process_post(i);
563 }
564 # endif
565
567 # endif
void ForSparseHesCheckpoint(vector< bool > &r, vector< bool > &s, local::sparse_list &h)
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
Add one element to a set.
size_t end(void) const
Fetch end for this vector of sets object.
File used to define pod_vector class.
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
void post_element(size_t i, size_t element)
Post an element for delayed addition to a set.
void ForSparseHesCase(bool set_type, const VectorSet &r, const VectorSet &s, VectorSet &h)
Private helper function for ForSparseHes(q, s) bool sparsity.
Two constant standard sets (currently used for concept checking).
VectorSet ForSparseHes(const VectorSet &r, const VectorSet &s)
User API for Hessian sparsity patterns using reverse mode.
Check that exp is true, if not terminate execution.
Add one element to a set.
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.
cons_iterator for one set of positive integers in a sparse_list object.
void rev_jac_sweep(const local::player< Base > *play, bool dependency, size_t n, size_t numvar, Vector_set &var_sparsity)
Given the sparsity pattern for the dependent variables, RevJacSweep computes the sparsity pattern for...
void for_hes_sweep(const local::player< Base > *play, size_t n, size_t numvar, const Vector_set &for_jac_sparse, const Vector_set &rev_jac_sparse, Vector_set &for_hes_sparse)
Given the forward Jacobian sparsity pattern for all the variables, and the reverse Jacobian sparsity ...
Scalar value_type
void for_jac_sweep(const local::player< Base > *play, bool dependency, size_t n, size_t numvar, Vector_set &var_sparsity)
Given the sparsity pattern for the independent variables, ForJacSweep computes the sparsity pattern f...
void resize(size_t n_set, size_t end)
Change number of sets, set end, and initialize all sets as empty.