CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
sparse_hes.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_CORE_SPARSE_HES_HPP
2 # define CPPAD_CORE_SPARSE_HES_HPP
3 /* --------------------------------------------------------------------------
4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
5 
6 CppAD is distributed under multiple licenses. This distribution is under
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.
11 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
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%
65  ADFun<%Base%> %f%
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 */
231 # include <cppad/core/cppad_assert.hpp>
235 
236 /*!
237 \file sparse_hes.hpp
238 Sparse Hessian calculation routines.
239 */
240 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
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
284 (this ADFun object).
285 
286 \param w
287 a vector of length m, the number of dependent variables in f
288 (this ADFun object).
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.
303 This must be cppad.symmetric, cppad.general, colpack.symmetic,
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
407  // (SizeVector and CppAD::vector<size_t>)
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" )
430  local::color_general_cppad(internal_pattern, col, row, color);
431  else if( coloring == "cppad.symmetric" )
432  local::color_symmetric_cppad(internal_pattern, col, row, color);
433  else if( coloring == "colpack.general" )
434  {
435 # if CPPAD_HAS_COLPACK
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  {
450 # if CPPAD_HAS_COLPACK
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  }
460  else CPPAD_ASSERT_KNOWN(
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 
538 } // END_CPPAD_NAMESPACE
539 
540 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
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.
CppAD::vector< size_t > order
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
CppAD::vector< size_t > col
Definition: sparse_hes.hpp:251
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
CppAD::vector< size_t > color
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
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
void color_symmetric_cppad(const VectorSet &pattern, CppAD::vector< size_t > &row, CppAD::vector< size_t > &col, CppAD::vector< size_t > &color)
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.
void color_symmetric_colpack(const VectorSet &pattern, CppAD::vector< size_t > &row, CppAD::vector< size_t > &col, CppAD::vector< size_t > &color)
Colpack algorithm for determining which rows of a symmetric sparse matrix can be computed together...
size_t nc(void) const
number of columns in matrix
Definition: sparse_rc.hpp:287
void set_internal_sparsity(bool zero_empty, bool input_empty, bool transpose, const vector< size_t > &internal_index, InternalSparsity &internal_pattern, const sparse_rc< SizeVector > &pattern_in)
Update the internal sparsity pattern for a sub-set of rows.
Routines that enable code to be independent of which internal spasity pattern is used.
const SizeVector & row(void) const
row indices
Definition: sparse_rcv.hpp:257
CppAD::vector< size_t > row
row and column indicies for return values (some may be reflected by symmetric coloring algorithms) ...
Definition: sparse_hes.hpp:250