\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

sparse_hes.cpp

View page source

Computing Sparse Hessian: Example and Test

# include <cppad/cppad.hpp>
bool sparse_hes(void)
{  bool ok = true;
   using CppAD::AD;
   using CppAD::NearEqual;
   //
   typedef CPPAD_TESTVECTOR(AD<double>)               a_vector;
   typedef CPPAD_TESTVECTOR(double)                   d_vector;
   typedef CPPAD_TESTVECTOR(size_t)                   s_vector;
   typedef CPPAD_TESTVECTOR(bool)                     b_vector;
   //
   // domain space vector
   size_t n = 12;  // must be greater than or equal 3; see n_sweep below
   a_vector a_x(n);
   for(size_t j = 0; j < n; j++)
      a_x[j] = AD<double> (0);
   //
   // declare independent variables and starting recording
   CppAD::Independent(a_x);
   //
   // range space vector
   size_t m = 1;
   a_vector a_y(m);
   a_y[0] = a_x[0] * a_x[1];
   for(size_t j = 0; j < n; j++)
      a_y[0] += a_x[j] * a_x[j] * a_x[j];
   //
   // create f: x -> y and stop tape recording
   // (without executing zero order forward calculation)
   CppAD::ADFun<double> f;
   f.Dependent(a_x, a_y);
   //
   // new value for the independent variable vector, and weighting vector
   d_vector w(m), x(n);
   for(size_t j = 0; j < n; j++)
      x[j] = double(j);
   w[0] = 1.0;
   //
   // vector used to check the value of the hessian
   d_vector check(n * n);
   size_t ij  = 0 * n + 1;
   for(ij = 0; ij < n * n; ij++)
      check[ij] = 0.0;
   ij         = 0 * n + 1;
   check[ij]  = 1.0;
   ij         = 1 * n + 0;
   check[ij]  = 1.0 ;
   for(size_t j = 0; j < n; j++)
   {  ij = j * n + j;
      check[ij] = 6.0 * x[j];
   }
   //
   // compute Hessian sparsity pattern
   b_vector select_domain(n), select_range(m);
   for(size_t j = 0; j < n; j++)
      select_domain[j] = true;
   select_range[0] = true;
   //
   CppAD::sparse_rc<s_vector> hes_pattern;
   bool internal_bool = false;
   f.for_hes_sparsity(
      select_domain, select_range, internal_bool, hes_pattern
   );
   //
   // compute entire sparse Hessian (really only need lower triangle)
   CppAD::sparse_rcv<s_vector, d_vector> subset( hes_pattern );
   CppAD::sparse_hes_work work;
   std::string coloring = "cppad.symmetric";
   size_t n_sweep = f.sparse_hes(x, w, subset, hes_pattern, coloring, work);
   ok &= n_sweep == 2;
   //
   const s_vector row( subset.row() );
   const s_vector col( subset.col() );
   const d_vector val( subset.val() );
   size_t nnz = subset.nnz();
   ok &= nnz == n + 2;
   for(size_t k = 0; k < nnz; k++)
   {  ij = row[k] * n + col[k];
      ok &= val[k] == check[ij];
   }
   return ok;
}