$\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}} }$
 # include <cppad/cppad.hpp> bool colpack_hessian(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::vector<size_t> i_vector; size_t i, j, k, ell; double eps = 10. * CppAD::numeric_limits<double>::epsilon(); // domain space vector size_t n = 5; 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); // colpack example case where hessian is a spear head // i.e, H(i, j) non zero implies i = 0, j = 0, or i = j AD<double> sum = 0.0; // partial_0 partial_j = x[j] // partial_j partial_j = x[0] for(j = 1; j < n; j++) sum += a_x[0] * a_x[j] * a_x[j] / 2.0; // // partial_i partial_i = 2 * x[i] for(i = 0; i < n; i++) sum += a_x[i] * a_x[i] * a_x[i] / 3.0; // declare dependent variables size_t m = 1; a_vector a_y(m); a_y[0] = sum; // create f: x -> y and stop tape recording CppAD::ADFun<double> f(a_x, a_y); // new value for the independent variable vector d_vector x(n); for(j = 0; j < n; j++) x[j] = double(j + 1); /* [ 2 2 3 4 5 ] hes = [ 2 5 0 0 0 ] [ 3 0 7 0 0 ] [ 4 0 0 9 0 ] [ 5 0 0 0 11 ] */ d_vector check(n * n); for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { size_t index = i * n + j; check[index] = 0.0; if( i == 0 && 1 <= j ) check[index] += x[j]; if( 1 <= i && j == 0 ) check[index] += x[i]; if( i == j ) { check[index] += 2.0 * x[i]; if( i != 0 ) check[index] += x[0]; } } } // Normally one would use f.RevSparseHes to compute // sparsity pattern, but for this example we extract it from check. std::vector< std::set<size_t> > p(n); i_vector row, col; for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { ell = i * n + j; if( check[ell] != 0. ) { // insert this non-zero entry in sparsity pattern p[i].insert(j); // the Hessian is symmetric, so only lower triangle if( j <= i ) { row.push_back(i); col.push_back(j); } } } } size_t K = row.size(); d_vector hes(K); // default coloring method is cppad.symmetric CppAD::sparse_hessian_work work; ok &= work.color_method == "cppad.symmetric"; // contrast and check results for both CppAD and Colpack for(size_t i_method = 0; i_method < 4; i_method++) { // empty work structure switch(i_method) { case 0: work.color_method = "cppad.symmetric"; break; case 1: work.color_method = "cppad.general"; break; case 2: work.color_method = "colpack.symmetric"; break; case 3: work.color_method = "colpack.general"; break; } // compute Hessian d_vector w(m); w[0] = 1.0; size_t n_sweep = f.SparseHessian(x, w, p, row, col, hes, work); // // check result for(k = 0; k < K; k++) { ell = row[k] * n + col[k]; ok &= NearEqual(check[ell], hes[k], eps, eps); } if( work.color_method == "cppad.symmetric" || work.color_method == "colpack.symmetric" ) ok &= n_sweep == 2; else ok &= n_sweep == 5; // // check that clear resets color_method to cppad.symmetric work.clear(); ok &= work.color_method == "cppad.symmetric"; } return ok; }