CppAD: A C++ Algorithmic Differentiation Package  20171217
sparse_hes.hpp
Go to the documentation of this file.
3 /* --------------------------------------------------------------------------
5
7 the terms of the
8  Eclipse Public License Version 1.0.
9
10 A copy of this license is included in the COPYING file of this distribution.
12 -------------------------------------------------------------------------- */
13
14 /*
15 $begin sparse_hes$$16 spell 17 const 18 Taylor 19 rc 20 rcv 21 nr 22 nc 23 hes 24 std 25 cppad 26 colpack 27 cmake 28 Jacobian 29$$ 30 31$section Computing Sparse Hessians$$32 33 head Syntax$$
34 $icode%n_sweep% = %f%.sparse_hes( 35 %x%, %w%, %subset%, %pattern%, %coloring%, %work% 36 )%$$37 38 head Purpose$$ 39 We use$latex F : \B{R}^n \rightarrow \B{R}^m$$to denote the 40 function corresponding to icode f$$.
41 Here $icode n$$is the cref/domain/seq_property/Domain/$$ size, 42 and$icode m$$is the cref/range/seq_property/Range/$$ size, or $icode f$$. 43 The syntax above takes advantage of sparsity when computing the Hessian 44 latex $45 H(x) = \dpow{2}{x} \sum_{i=0}^{m-1} w_i F_i (x) 46$$$ 47 In the sparse case, this should be faster and take less memory than 48$cref Hessian$$. 49 The matrix element latex H_{i,j} (x)$$ is the second partial of
50 $latex w^\R{T} F (x)$$with respect to latex x_i$$ and$latex x_j$$. 51 52 head SizeVector$$
53 The type $icode SizeVector$$is a cref SimpleVector$$ class with 54$cref/elements of type/SimpleVector/Elements of Specified Type/$$55 code size_t$$.
56
57 $head BaseVector$$58 The type icode BaseVector$$ 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 This object has prototype 64$codei%
66 %$$67 Note that the Taylor coefficients stored in icode f$$ are affected
68 by this operation; see
69 $cref/uses forward/sparse_hes/Uses Forward/$$below. 70 71 head x$$ 72 This argument has prototype 73$codei%
74  const %BaseVector%& %x%
75 %$$76 and its size is icode n$$.
77 It specifies the point at which to evaluate the Hessian
78 $latex H(x)$$. 79 80 head w$$ 81 This argument has prototype 82$codei%
83  const %BaseVector%& %w%
84 %$$85 and its size is icode m$$.
86 It specifies the weight for each of the components of $latex F(x)$$; 87 i.e. latex w_i$$ is the weight for$latex F_i (x)$$. 88 89 head subset$$
90 This argument has prototype
91 $codei% 92 sparse_rcv<%SizeVector%, %BaseVector%>& %subset% 93 %$$94 Its row size and column size is icode n$$; i.e., 95$icode%subset%.nr() == %n%$$and icode%subset%.nc() == %n%$$.
96 It specifies which elements of the Hessian are computed.
97 $list number$$98 The input value of its value vector 99 icode%subset%.val()%$$ does not matter. 100 Upon return it contains the value of the corresponding elements 101 of the Hessian. 102$lnext
103 All of the row, column pairs in $icode subset$$must also appear in 104 icode pattern$$; i.e., they must be possibly non-zero. 105$lnext
106 The Hessian is symmetric, so one has a choice as to which off diagonal
107 elements to put in $icode subset$$. 108 It will probably be more efficient if one makes this choice so that 109 the there are more entries in each non-zero column of icode subset$$; 110 see$cref/n_sweep/sparse_hes/n_sweep/$$below. 111 lend 112 113 head pattern$$
114 This argument has prototype
115 $codei% 116 const sparse_rc<%SizeVector%>& %pattern% 117 %$$118 Its row size and column size is icode n$$; i.e., 119$icode%pattern%.nr() == %n%$$and icode%pattern%.nc() == %n%$$.
120 It is a sparsity pattern for the Hessian $latex H(x)$$. 121 If the th i$$ row ($th j$$column) does not appear in icode subset$$,
122 the $th i$$row (th j$$ column) of$icode pattern$$does not matter 123 and need not be computed. 124 This argument is not used (and need not satisfy any conditions), 125 when cref/work/sparse_hes/work/$$ is non-empty.
126
127 $subhead subset$$128 If the th i$$ row and$th i$$column do not appear in icode subset$$,
129 the $th i$$row and column of icode pattern$$ do not matter. 130 In this case the$th i-th$$row and column may have no entries in 131 icode pattern$$ even though they are possibly non-zero in $latex H(x)$$. 132 (This can be used to reduce the amount of computation required to find 133 icode pattern$$.) 134 135$head coloring$$136 The coloring algorithm determines which rows and columns 137 can be computed during the same sweep. 138 This field has prototype 139 codei% 140 const std::string& %coloring% 141 %$$
142 This value only matters when work is empty; i.e.,
143 after the $icode work$$constructor or icode%work%.clear()%$$. 144 145$subhead cppad.symmetric$$146 This coloring takes advantage of the fact that the Hessian matrix 147 is symmetric when find a coloring that requires fewer 148 cref/sweeps/sparse_hes/n_sweep/$$.
149
150 $subhead cppad.general$$151 This is the same as the sparse Jacobian 152 cref/cppad/sparse_jac/coloring/cppad/$$ method 153 which does not take advantage of symmetry. 154 155$subhead colpack.symmetric$$156 If cref colpack_prefix$$ was specified on the
157 $cref/cmake command/cmake/CMake Command/$$line, 158 you can set icode coloring$$ to$code colpack.symmetric$$. 159 This also takes advantage of the fact that the Hessian matrix is symmetric. 160 161 subhead colpack.general$$
162 If $cref colpack_prefix$$was specified on the 163 cref/cmake command/cmake/CMake Command/$$ line, 164 you can set$icode coloring$$to code colpack.general$$.
165 This is the same as the sparse Jacobian
166 $cref/colpack/sparse_jac/coloring/colpack/$$method 167 which does not take advantage of symmetry. 168 169 subhead colpack.star Deprecated 2017-06-01$$ 170 The$code colpack.star$$method is deprecated. 171 It is the same as the code colpack.symmetric$$ method
172 which should be used instead.
173
174
175 $head work$$176 This argument has prototype 177 codei% 178 sparse_hes_work& %work% 179 %$$ 180 We refer to its initial value, 181 and its value after$icode%work%.clear()%$$, as empty. 182 If it is empty, information is stored in icode work$$.
183 This can be used to reduce computation when
184 a future call is for the same object $icode f$$, 185 and the same subset of the Hessian. 186 If either of these values change, use icode%work%.clear()%$$ to 187 empty this structure. 188 189$head n_sweep$$190 The return value icode n_sweep$$ has prototype
191 $codei% 192 size_t %n_sweep% 193 %$$194 It is the number of first order forward sweeps 195 used to compute the requested Hessian values. 196 Each first forward sweep is followed by a second order reverse sweep 197 so it is also the number of reverse sweeps. 198 It is also the number of colors determined by the coloring method 199 mentioned above. 200 This is proportional to the total computational work, 201 not counting the zero order forward sweep, 202 or combining multiple columns and rows into a single sweep. 203 204 head Uses Forward$$ 205 After each call to$cref Forward$$, 206 the object icode f$$ contains the corresponding
207 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 208 After a call to code sparse_hes$$ 209 the zero order coefficients correspond to 210$codei%
211  %f%.Forward(0, %x%)
212 %$$213 All the other forward mode coefficients are unspecified. 214 215 head Example$$
216 $children% 217 example/sparse/sparse_hes.cpp 218 %$$219 The files cref sparse_hes.cpp$$ 220 is an example and test of$code sparse_hes$$. 221 It returns code true$$, if it succeeds, and $code false$$otherwise. 222 223 head Subset Hessian$$ 224 The routine 225$cref sparse_sub_hes.cpp$$226 is an example and test that compute a subset of a sparse Hessian. 227 It returns code true$$, for success, and $code false$$otherwise. 228 229$end
230 */
235
236 /*!
237 \file sparse_hes.hpp
238 Sparse Hessian calculation routines.
239 */
241
242 /*!
243 Class used to hold information used by Sparse Hessian routine in this file,
244 so it does not need to be recomputed every time.
245 */
247  public:
248  /// row and column indicies for return values
249  /// (some may be reflected by symmetric coloring algorithms)
252  /// indices that sort the row and col arrays by color
254  /// results of the coloring algorithm
256
257  /// constructor
259  { }
260  /// inform CppAD that this information needs to be recomputed
261  void clear(void)
262  {
263  row.clear();
264  col.clear();
265  order.clear();
266  color.clear();
267  }
268 };
269 // ----------------------------------------------------------------------------
270 /*!
271 Calculate sparse Hessians using forward mode
272
273 \tparam Base
274 the base type for the recording that is stored in the ADFun object.
275
276 \tparam SizeVector
277 a simple vector class with elements of type size_t.
278
279 \tparam BaseVector
280 a simple vector class with elements of type Base.
281
282 \param x
283 a vector of length n, the number of independent variables in f
285
286 \param w
287 a vector of length m, the number of dependent variables in f
289
290 \param subset
291 specifices the subset of the sparsity pattern where the Hessian is evaluated.
292 subset.nr() == n,
293 subset.nc() == n.
294
295 \param pattern
296 is a sparsity pattern for the Hessian of w^T * f;
297 pattern.nr() == n,
298 pattern.nc() == n,
299 where m is number of dependent variables in f.
300
301 \param coloring
302 determines which coloring algorithm is used.
304 or colpack.star.
305
306 \param work
307 this structure must be empty, or contain the information stored
308 by a previous call to sparse_hes.
309 The previous call must be for the same ADFun object f
310 and the same subset.
311
312 \return
313 This is the number of first order forward
314 (and second order reverse) sweeps used to compute thhe Hessian.
315 */
316 template <class Base>
317 template <class SizeVector, class BaseVector>
319  const BaseVector& x ,
320  const BaseVector& w ,
322  const sparse_rc<SizeVector>& pattern ,
323  const std::string& coloring ,
324  sparse_hes_work& work )
325 { size_t n = Domain();
326  //
328  subset.nr() == n,
329  "sparse_hes: subset.nr() not equal domain dimension for f"
330  );
332  subset.nc() == n,
333  "sparse_hes: subset.nc() not equal domain dimension for f"
334  );
336  size_t( x.size() ) == n,
337  "sparse_hes: x.size() not equal domain dimension for f"
338  );
340  size_t( w.size() ) == Range(),
341  "sparse_hes: w.size() not equal range dimension for f"
342  );
343  //
344  // work information
345  vector<size_t>& row(work.row);
346  vector<size_t>& col(work.col);
347  vector<size_t>& color(work.color);
348  vector<size_t>& order(work.order);
349  //
350  // subset information
351  const SizeVector& subset_row( subset.row() );
352  const SizeVector& subset_col( subset.col() );
353  //
354  // point at which we are evaluationg the Hessian
355  Forward(0, x);
356  //
357  // number of elements in the subset
358  size_t K = subset.nnz();
359  //
360  // check for case were there is nothing to do
361  // (except for call to Forward(0, x)
362  if( K == 0 )
363  return 0;
364  //
365 # ifndef NDEBUG
366  if( color.size() != 0 )
368  color.size() == n,
369  "sparse_hes: work is non-empty and conditions have changed"
370  );
372  row.size() == K,
373  "sparse_hes: work is non-empty and conditions have changed"
374  );
376  col.size() == K,
377  "sparse_hes: work is non-empty and conditions have changed"
378  );
379  //
380  for(size_t k = 0; k < K; k++)
381  { bool ok = row[k] == subset_row[k] && col[k] == subset_col[k];
382  ok |= row[k] == subset_col[k] && col[k] == subset_row[k];
384  ok,
385  "sparse_hes: work is non-empty and conditions have changed"
386  );
387  }
388  }
389 # endif
390  //
391  // check for case where input work is empty
392  if( color.size() == 0 )
393  { // compute work color and order vectors
395  pattern.nr() == n,
396  "sparse_hes: pattern.nr() not equal domain dimension for f"
397  );
399  pattern.nc() == n,
400  "sparse_hes: pattern.nc() not equal domain dimension for f"
401  );
402  //
403  // initialize work row, col to be same as subset row, col
404  row.resize(K);
405  col.resize(K);
406  // cannot assign vectors becasue may be of different types
408  for(size_t k = 0; k < K; k++)
409  { row[k] = subset_row[k];
410  col[k] = subset_col[k];
411  }
412  //
413  // convert pattern to an internal version of its transpose
414  vector<size_t> internal_index(n);
415  for(size_t j = 0; j < n; j++)
416  internal_index[j] = j;
417  bool transpose = true;
418  bool zero_empty = false;
419  bool input_empty = true;
420  local::sparse_list internal_pattern;
421  internal_pattern.resize(n, n);
422  local::set_internal_sparsity(zero_empty, input_empty,
423  transpose, internal_index, internal_pattern, pattern
424  );
425  //
426  // execute coloring algorithm
427  // (we are using transpose becasue coloring groups rows, not columns)
428  color.resize(n);
429  if( coloring == "cppad.general" )
431  else if( coloring == "cppad.symmetric" )
433  else if( coloring == "colpack.general" )
434  {
436  local::color_general_colpack(internal_pattern, col, row, color);
437 # else
439  false,
440  "sparse_hes: coloring = colpack.star "
441  "and colpack_prefix not in cmake command line."
442  );
443 # endif
444  }
445  else if(
446  coloring == "colpack.symmetric" ||
447  coloring == "colpack.star"
448  )
449  {
451  local::color_symmetric_colpack(internal_pattern, col, row, color);
452 # else
454  false,
455  "sparse_hes: coloring = colpack.symmetic or colpack.star "
456  "and colpack_prefix not in cmake command line."
457  );
458 # endif
459  }
461  false,
462  "sparse_hes: coloring is not valid."
463  );
464  //
465  // put sorting indices in color order
466  SizeVector key(K);
467  order.resize(K);
468  for(size_t k = 0; k < K; k++)
469  key[k] = color[ col[k] ];
470  index_sort(key, order);
471  }
472  // Base versions of zero and one
473  Base one(1.0);
474  Base zero(0.0);
475  //
476  size_t n_color = 1;
477  for(size_t j = 0; j < n; j++) if( color[j] < n )
478  n_color = std::max(n_color, color[j] + 1);
479  //
480  // initialize the return Hessian values as zero
481  for(size_t k = 0; k < K; k++)
482  subset.set(k, zero);
483  //
484  // direction vector for calls to first order forward
485  BaseVector dx(n);
486  //
487  // return values for calls to second order reverse
488  BaseVector ddw(2 * n);
489  //
490  // loop over colors
491  size_t k = 0;
492  for(size_t ell = 0; ell < n_color; ell++)
493  if( k == K )
494  { // kludge because colpack returns colors that are not used
495  // (it does not know about the subset corresponding to row, col)
497  coloring == "colpack.general" ||
498  coloring == "colpack.symmetric" ||
499  coloring == "colpack.star"
500  );
501  }
502  else if( color[ col[ order[k] ] ] != ell )
503  { // kludge because colpack returns colors that are not used
504  // (it does not know about the subset corresponding to row, col)
506  coloring == "colpack.general" ||
507  coloring == "colpack.symmetic" ||
508  coloring == "colpack.star"
509  );
510  }
511  else
512  { CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
513  //
514  // combine all columns with this color
515  for(size_t j = 0; j < n; j++)
516  { dx[j] = zero;
517  if( color[j] == ell )
518  dx[j] = one;
519  }
520  // call forward mode for all these rows at once
521  Forward(1, dx);
522  //
523  // evaluate derivative of w^T * F'(x) * dx
524  ddw = Reverse(2, w);
525  //
526  // set the corresponding components of the result
527  while( k < K && color[ col[order[k]] ] == ell )
528  { size_t index = row[ order[k] ] * 2 + 1;
529  subset.set(order[k], ddw[index] );
530  k++;
531  }
532  }
533  // check that all the required entries have been set
534  CPPAD_ASSERT_UNKNOWN( k == K );
535  return n_color;
536 }
537
539
540 # 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
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 set(size_t k, const value_type &v)
Definition: sparse_rcv.hpp:240
size_t nnz(void) const
number of possibly non-zero elements in matrix
Definition: sparse_rcv.hpp:254
const SizeVector & col(void) const
column indices
Definition: sparse_rcv.hpp:260
sparse_hes_work(void)
constructor
Definition: sparse_hes.hpp:258
Coloring algorithm for a symmetric sparse matrix.
indices that sort the row and col arrays by color
Definition: sparse_hes.hpp:253
void clear(void)
free memory and set number of elements to zero
Definition: vector.hpp:418
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
size_t nc(void) const
number of columns in matrix
Definition: sparse_rcv.hpp:251
Sparse matrices with elements of type Scalar.
Definition: sparse_rcv.hpp:212
Coloring algorithm for a general sparse matrix.
void resize(size_t n)
change the number of elements in this vector.
Definition: vector.hpp:399
size_t nr(void) const
number of rows in matrix
Definition: sparse_rcv.hpp:248
size_t nr(void) const
number of rows in matrix
Definition: sparse_rc.hpp:284
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
Definition: sparse_hes.hpp:251
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
results of the coloring algorithm
Definition: sparse_hes.hpp:255
void index_sort(const VectorKey &keys, VectorSize &ind)
Compute the indices that sort a vector of keys.
Definition: index_sort.hpp:140
Class used to hold information used by Sparse Hessian routine in this file, so it does not need to be...
Definition: sparse_hes.hpp:246
void clear(void)
inform CppAD that this information needs to be recomputed
Definition: sparse_hes.hpp:261
size_t sparse_hes(const BaseVector &x, const BaseVector &w, sparse_rcv< SizeVector, BaseVector > &subset, const sparse_rc< SizeVector > &pattern, const std::string &coloring, sparse_hes_work &work)
Calculate sparse Hessians using forward mode.
Definition: sparse_hes.hpp:318
Check that exp is true, if not terminate execution.
CppAD algorithm for determining which rows of a symmetric sparse matrix can be computed together...
void resize(size_t n_set, size_t end)
Start a new vector of sets.