CppAD: A C++ Algorithmic Differentiation Package  20171217
for_hes_sparsity.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 $begin for_hes_sparsity$$16 spell 17 Andrea Walther 18 Jacobian 19 Hessian 20 jac 21 hes 22 bool 23 const 24 rc 25 cpp 26$$ 27 28$section Forward Mode Hessian Sparsity Patterns$$29 30 head Syntax$$
31 $icode%f%.for_hes_sparsity( 32 %select_domain%, %select_range%, %internal_bool%, %pattern_out% 33 )%$$34 35 head Purpose$$ 36 We use$latex F : \B{R}^n \rightarrow \B{R}^m$$to denote the 37 cref/AD function/glossary/AD Function/$$ corresponding to
38 the operation sequence stored in $icode f$$. 39 Fix a diagonal matrix latex D \in \B{R}^{n \times n}$$, 40 a vector$latex s \in \B{R}^m$$and define the function 41 latex $42 H(x) = D ( s^\R{T} F )^{(2)} ( x ) D 43$$$
44 Given the sparsity for $latex D$$and latex s$$, 45$code for_hes_sparsity$$computes a sparsity pattern for latex H(x)$$.
46
47 $head x$$48 Note that the sparsity pattern latex H(x)$$ corresponds to the 49 operation sequence stored in$icode f$$and does not depend on 50 the argument icode x$$.
51
52 $head BoolVector$$53 The type icode BoolVector$$ is a$cref SimpleVector$$class with 54 cref/elements of type/SimpleVector/Elements of Specified Type/$$
55 $code bool$$. 56 57 head SizeVector$$ 58 The type$icode SizeVector$$is a cref SimpleVector$$ class with
59 $cref/elements of type/SimpleVector/Elements of Specified Type/$$60 code size_t$$. 61 62$head f$$63 The object icode f$$ has prototype
64 $codei% 65 ADFun<%Base%> %f% 66 %$$67 68 head select_domain$$ 69 The argument$icode select_domain$$has prototype 70 codei% 71 const %BoolVector%& %select_domain% 72 %$$
73 It has size $latex n$$and specifies which components of the diagonal of 74 latex D$$ are non-zero; i.e.,$icode%select_domain%[%j%]%$$is true 75 if and only if latex D_{j,j}$$ is possibly non-zero.
76
77
78 $head select_range$$79 The argument icode select_range$$ has prototype 80$codei%
81  const %BoolVector%& %select_range%
82 %$$83 It has size latex m$$ and specifies which components of the vector
84 $latex s$$are non-zero; i.e., icode%select_range%[%i%]%$$ is true 85 if and only if$latex s_i$$is possibly non-zero. 86 87 head internal_bool$$
88 If this is true, calculations are done with sets represented by a vector
89 of boolean values. Otherwise, a vector of sets of integers is used.
90
91 $head pattern_out$$92 This argument has prototype 93 codei% 94 sparse_rc<%SizeVector%>& %pattern_out% 95 %$$ 96 This input value of$icode pattern_out$$does not matter. 97 Upon return icode pattern_out$$ is a sparsity pattern for $latex H(x)$$. 98 99 head Sparsity for Entire Hessian$$ 100 Suppose that$latex R$$is the latex n \times n$$ identity matrix.
101 In this case, $icode pattern_out$$is a sparsity pattern for 102 latex (s^\R{T} F) F^{(2)} ( x )$$. 103 104$head Algorithm$$105 See Algorithm II in 106 italic Computing sparse Hessians with automatic differentiation$$
107 by Andrea Walther.
108 Note that $icode s$$provides the information so that 109 'dead ends' are not included in the sparsity pattern. 110 111 head Example$$ 112$children%
113  example/sparse/for_hes_sparsity.cpp
114 %$$115 The file cref for_hes_sparsity.cpp$$
116 contains an example and test of this operation.
117 It returns true if it succeeds and false otherwise.
118
119 \$end
120 -----------------------------------------------------------------------------
121 */
124
126
127 /*!
128 Forward Hessian sparsity patterns.
129
130 \tparam Base
131 is the base type for this recording.
132
133 \tparam BoolVector
134 is the simple vector with elements of type bool that is used for
135 sparsity for the vector s.
136
137 \tparam SizeVector
138 is the simple vector with elements of type size_t that is used for
139 row, column index sparsity patterns.
140
141 \param select_domain
142 is a sparsity pattern for for the diagonal of D.
143
144 \param select_range
145 is a sparsity pattern for for s.
146
147 \param internal_bool
148 If this is true, calculations are done with sets represented by a vector
149 of boolean values. Otherwise, a vector of standard sets is used.
150
151 \param pattern_out
152 The return value is a sparsity pattern for H(x) where
153 \f[
154  H(x) = D * F^{(1)} (x) * D
155 \f]
156 Here F is the function corresponding to the operation sequence
157 and x is any argument value.
158 */
159 template <class Base>
160 template <class BoolVector, class SizeVector>
162  const BoolVector& select_domain ,
163  const BoolVector& select_range ,
164  bool internal_bool ,
165  sparse_rc<SizeVector>& pattern_out )
166 { size_t n = Domain();
167  size_t m = Range();
168  //
170  size_t( select_domain.size() ) == n,
171  "for_hes_sparsity: size of select_domain is not equal to "
172  "number of independent variables"
173  );
175  size_t( select_range.size() ) == m,
176  "for_hes_sparsity: size of select_range is not equal to "
177  "number of dependent variables"
178  );
179  // do not need transpose or depenency
180  bool transpose = false;
181  bool dependency = false;
182  //
183  sparse_rc<SizeVector> pattern_tmp;
184  if( internal_bool )
185  { // forward Jacobian sparsity pattern for independent variables
186  local::sparse_pack internal_for_jac;
187  internal_for_jac.resize(num_var_tape_, n + 1 );
188  for(size_t j = 0; j < n; j++) if( select_domain[j] )
192  }
193  // forward Jacobian sparsity for all variables on tape
195  &play_,
196  dependency,
197  n,
198  num_var_tape_,
199  internal_for_jac
200  );
201  // reverse Jacobian sparsity pattern for select_range
202  local::sparse_pack internal_rev_jac;
203  internal_rev_jac.resize(num_var_tape_, 1);
204  for(size_t i = 0; i < m; i++) if( select_range[i] )
208  }
209  // reverse Jacobian sparsity for all variables on tape
211  &play_,
212  dependency,
213  n,
214  num_var_tape_,
215  internal_rev_jac
216  );
217  // internal vector of sets that will hold Hessian
218  local::sparse_pack internal_for_hes;
219  internal_for_hes.resize(n + 1, n + 1);
220  //
221  // compute forward Hessian sparsity pattern
223  &play_,
224  n,
225  num_var_tape_,
226  internal_for_jac,
227  internal_rev_jac,
228  internal_for_hes
229  );
230  //
231  // put the result in pattern_tmp
234  );
235  }
236  else
237  { // forward Jacobian sparsity pattern for independent variables
238  // (corresponds to D)
239  local::sparse_list internal_for_jac;
240  internal_for_jac.resize(num_var_tape_, n + 1 );
241  for(size_t j = 0; j < n; j++) if( select_domain[j] )
245  }
246  // forward Jacobian sparsity for all variables on tape
248  &play_,
249  dependency,
250  n,
251  num_var_tape_,
252  internal_for_jac
253  );
254  // reverse Jacobian sparsity pattern for select_range
255  // (corresponds to s)
256  local::sparse_list internal_rev_jac;
257  internal_rev_jac.resize(num_var_tape_, 1);
258  for(size_t i = 0; i < m; i++) if( select_range[i] )
262  }
263  // reverse Jacobian sparsity for all variables on tape
265  &play_,
266  dependency,
267  n,
268  num_var_tape_,
269  internal_rev_jac
270  );
271  // internal vector of sets that will hold Hessian
272  local::sparse_list internal_for_hes;
273  internal_for_hes.resize(n + 1, n + 1);
274  //
275  // compute forward Hessian sparsity pattern
277  &play_,
278  n,
279  num_var_tape_,
280  internal_for_jac,
281  internal_rev_jac,
282  internal_for_hes
283  );
284  //
285  // put the result in pattern_tmp
288  );
289  }
290  // subtract 1 from all column values
291  CPPAD_ASSERT_UNKNOWN( pattern_tmp.nr() == n );
292  CPPAD_ASSERT_UNKNOWN( pattern_tmp.nc() == n + 1 );
293  const SizeVector& row( pattern_tmp.row() );
294  const SizeVector& col( pattern_tmp.col() );
295  size_t nr = n;
296  size_t nc = n;
297  size_t nnz = pattern_tmp.nnz();
298  pattern_out.resize(nr, nc, nnz);
299  for(size_t k = 0; k < nnz; k++)
300  { CPPAD_ASSERT_UNKNOWN( 0 < col[k] );
301  pattern_out.set(k, row[k], col[k] - 1);
302  }
303  return;
304 }
306 # endif
Check that exp is true, if not print msg and terminate execution.
Vector of sets of positive integers, each set stored as a singly linked list.
Definition: sparse_list.hpp:35
const SizeVector & col(void) const
column indices
Definition: sparse_rc.hpp:296
void for_hes_sparsity(const BoolVector &select_domain, const BoolVector &select_range, bool internal_bool, sparse_rc< SizeVector > &pattern_out)
Forward Hessian sparsity patterns.
Add one element to a set.
void get_internal_sparsity(bool transpose, const vector< size_t > &internal_index, const InternalSparsity &internal_pattern, sparse_rc< SizeVector > &pattern_out)
Get sparsity pattern for a sub-set of variables.
size_t nr(void) const
number of rows in matrix
Definition: sparse_rc.hpp:284
Vector of sets of postivie integers, each set stored as a packed boolean array.
Definition: sparse_pack.hpp:33
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
size_t nnz(void) const
number of possibly non-zero elements in matrix
Definition: sparse_rc.hpp:290
void resize(size_t nr, size_t nc, size_t nnz)
resize
Definition: sparse_rc.hpp:258
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.
size_t nc(void) const
number of columns in matrix
Definition: sparse_rc.hpp:287
File used to define the ADFun&lt;Base&gt; class.
void set(size_t k, size_t r, size_t c)
set row and column for a possibly non-zero element
Definition: sparse_rc.hpp:266
Routines that enable code to be independent of which internal spasity pattern is used.
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 ...
const SizeVector & row(void) const
row indices
Definition: sparse_rc.hpp:293
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.