$\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_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::vector<size_t> i_vector; typedef CppAD::sparse_rc<i_vector> sparsity; typedef CppAD::sparse_rcv<i_vector, d_vector> sparse_matrix; double eps = 10. * CppAD::numeric_limits<double>::epsilon(); // // domain space vector size_t n = 5; 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); // 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(size_t 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(size_t 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(size_t 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 ] */ // Normally one would use CppAD to compute sparsity pattern, but for this // example we set it directly size_t nr = n; size_t nc = n; size_t nnz = n + 2 * (n - 1); sparsity pattern(nr, nc, nnz); for(size_t k = 0; k < n; k++) { size_t r = k; size_t c = k; pattern.set(k, r, c); } for(size_t i = 1; i < n; i++) { size_t k = n + 2 * (i - 1); size_t r = i; size_t c = 0; pattern.set(k, r, c); pattern.set(k+1, c, r); } // subset of elements to compute // (only compute lower traingle) nnz = n + (n - 1); sparsity lower_triangle(nr, nc, nnz); d_vector check(nnz); for(size_t k = 0; k < n; k++) { size_t r = k; size_t c = k; lower_triangle.set(k, r, c); check[k] = 2.0 * x[k]; if( k > 0 ) check[k] += x[0]; } for(size_t j = 1; j < n; j++) { size_t k = n + (j - 1); size_t r = 0; size_t c = j; lower_triangle.set(k, r, c); check[k] = x[c]; } sparse_matrix subset( lower_triangle ); // check results for both CppAD and Colpack for(size_t i_method = 0; i_method < 4; i_method++) { // coloring method std::string coloring; switch(i_method) { case 0: coloring = "cppad.symmetric"; break; case 1: coloring = "cppad.general"; break; case 2: coloring = "colpack.symmetric"; break; case 3: coloring = "colpack.general"; break; } // // compute Hessian CppAD::sparse_hes_work work; d_vector w(m); w[0] = 1.0; size_t n_sweep = f.sparse_hes( x, w, subset, pattern, coloring, work ); // // check result const d_vector& hes( subset.val() ); for(size_t k = 0; k < nnz; k++) ok &= NearEqual(check[k], hes[k], eps, eps); if( coloring == "cppad.symmetric" || coloring == "colpack.symmetric" ) ok &= n_sweep == 2; else ok &= n_sweep == 5; } return ok; }