sparse_hessian.cpp¶

Sparse Hessian: Example and Test¶

# include <cppad/cppad.hpp>
bool sparse_hessian(void)
{  bool ok = true;
   using CppAD::AD;
   using CppAD::NearEqual;
   size_t i, j, k, ell;
   typedef CPPAD_TESTVECTOR(AD<double>)               a_vector;
   typedef CPPAD_TESTVECTOR(double)                     d_vector;
   typedef CPPAD_TESTVECTOR(size_t)                     i_vector;
   typedef CPPAD_TESTVECTOR(bool)                       b_vector;
   typedef CPPAD_TESTVECTOR(std::set<size_t>)         s_vector;
   double eps = 10. * CppAD::numeric_limits<double>::epsilon();

   // domain space vector
   size_t n = 12;  // must be greater than or equal 3; see n_sweep below
   a_vector a_x(n);
   for(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(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(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);
   for(ell = 0; ell < n * n; ell++)
      check[ell] = 0.0;
   ell        = 0 * n + 1;
   check[ell] = 1.0;
   ell        = 1 * n + 0;
   check[ell] = 1.0 ;
   for(j = 0; j < n; j++)
   {  ell = j * n + j;
      check[ell] = 6.0 * x[j];
   }

   // -------------------------------------------------------------------
   // second derivative of y[0] w.r.t x
   d_vector hes(n * n);
   hes = f.SparseHessian(x, w);
   for(ell = 0; ell < n * n; ell++)
      ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

   // --------------------------------------------------------------------
   // example using vectors of bools to compute sparsity pattern for Hessian
   b_vector r_bool(n * n);
   for(i = 0; i < n; i++)
   {  for(j = 0; j < n; j++)
         r_bool[i * n + j] = false;
      r_bool[i * n + i] = true;
   }
   f.ForSparseJac(n, r_bool);
   //
   b_vector s_bool(m);
   for(i = 0; i < m; i++)
      s_bool[i] = w[i] != 0;
   b_vector p_bool = f.RevSparseHes(n, s_bool);

   hes = f.SparseHessian(x, w, p_bool);
   for(ell = 0; ell < n * n; ell++)
      ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

   // --------------------------------------------------------------------
   // example using vectors of sets to compute sparsity pattern for Hessian
   s_vector r_set(n);
   for(i = 0; i < n; i++)
      r_set[i].insert(i);
   f.ForSparseJac(n, r_set);
   //
   s_vector s_set(m);
   for(i = 0; i < m; i++)
      if( w[i] != 0. )
         s_set[0].insert(i);
   s_vector p_set = f.RevSparseHes(n, s_set);

   // example passing sparsity pattern to SparseHessian
   hes = f.SparseHessian(x, w, p_set);
   for(ell = 0; ell < n * n; ell++)
      ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

   // --------------------------------------------------------------------
   // use row and column indices to specify upper triangle of
   // non-zero elements of Hessian
   size_t K = n + 1;
   i_vector row(K), col(K);
   hes.resize(K);
   k = 0;
   for(j = 0; j < n; j++)
   {  // diagonal of Hessian
      row[k] = j;
      col[k] = j;
      k++;
   }
   // only off diagonal non-zero elemenet in upper triangle
   row[k] = 0;
   col[k] = 1;
   k++;
   ok &= k == K;
   CppAD::sparse_hessian_work work;

   // can use p_set or p_bool.
   size_t n_sweep = f.SparseHessian(x, w, p_set, row, col, hes, work);
   for(k = 0; k < K; k++)
   {  ell = row[k] * n + col[k];
      ok &=  NearEqual(w[0] * check[ell], hes[k], eps, eps );
   }
   ok &= n_sweep == 2;

   // now recompute at a different x and w (using work from previous call
   w[0]       = 2.0;
   x[1]       = 0.5;
   ell        = 1 * n + 1;
   check[ell] = 6.0 * x[1];
   s_vector   not_used;
   n_sweep    = f.SparseHessian(x, w, not_used, row, col, hes, work);
   for(k = 0; k < K; k++)
   {  ell = row[k] * n + col[k];
      ok &=  NearEqual(w[0] * check[ell], hes[k], eps, eps );
   }
   ok &= n_sweep == 2;



   return ok;
}