$\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}} }$
CppAD Speed: Sparse Hessian

Specifications

Implementation
# include <cppad/cppad.hpp>

// Note that CppAD uses global_option["memory"] at the main program level
# include <map>
extern std::map<std::string, bool> global_option;

namespace {
// typedefs
typedef vector<bool>                          b_vector;
typedef vector<size_t>                        s_vector;
typedef vector<double>                        d_vector;
typedef vector<a1double>                      a1vector;
typedef vector<a2double>                      a2vector;
typedef CppAD::sparse_rcv<s_vector, d_vector> sparse_matrix;
// ------------------------------------------------------------------------
void create_fun(
const d_vector&             x        ,
const s_vector&             row      ,
const s_vector&             col      ,
{
// initialize a1double version of independent variables
size_t n = x.size();
a1vector a1x(n);
for(size_t j = 0; j < n; j++)
a1x[j] = x[j];
//
// optimization options
std::string optimize_options="no_compare_op";
//
// order of derivative in sparse_hes_fun
size_t order = 0;
//
if( ! global_option["hes2jac"] )
{
// declare independent variables
Independent(a1x);
//
// AD computation of y
a1vector a1y(1);
CppAD::sparse_hes_fun<a1double>(n, a1x, row, col, order, a1y);
//
// create function object f : X -> Y
fun.Dependent(a1x, a1y);
//
if( global_option["optimize"] )
fun.optimize(optimize_options);
//
// skip comparison operators
fun.compare_change_count(0);
//
// fun corresonds to f(x)
return;
}
// declare independent variables for f(x)
a2vector a2x(n);
for(size_t j = 0; j < n; j++)
a2x[j] = a1x[j];
Independent(a2x);
//
// a2double computation of y
a2vector a2y(1);
CppAD::sparse_hes_fun<a2double>(n, a2x, row, col, order, a2y);
//
// create function object corresponding to y = f(x)
a1f.Dependent(a2x, a2y);
//
// declare independent variables for g(x)
Independent(a1x);
//
// a1double computation of z
a1vector a1w(1), a1z(n);
a1w[0] = 1.0;
a1f.Forward(0, a1x);
a1z = a1f.Reverse(1, a1w);
//
// create function object z = g(x) = f'(x)
fun.Dependent(a1x, a1z);
//
if( global_option["optimize"] )
fun.optimize(optimize_options);
//
// skip comparison operators
fun.compare_change_count(0);
//
// fun corresonds to g(x)
return;
}
// ------------------------------------------------------------------------
void calc_sparsity(
sparsity_pattern&      sparsity ,
{
size_t n = fun.Domain();
size_t m = fun.Range();
//
bool transpose     = false;
//
if( global_option["subsparsity"] )
{     CPPAD_ASSERT_UNKNOWN( global_option["hes2jac"] )
CPPAD_ASSERT_UNKNOWN( n == m );
b_vector select_domain(n), select_range(m);
for(size_t j = 0; j < n; ++j)
select_domain[j] = true;
for(size_t i = 0; i < m; ++i)
select_range[i] = true;
//
// fun corresponds to g(x)
fun.subgraph_sparsity(
select_domain, select_range, transpose, sparsity
);
return;
}
bool dependency    = false;
bool reverse       = global_option["revsparsity"];
bool internal_bool = global_option["boolsparsity"];
//
if( ! global_option["hes2jac"] )
{     // fun corresponds to f(x)
//
CPPAD_ASSERT_UNKNOWN( m == 1 );
//
b_vector select_range(m);
select_range[0] = true;
//
if( reverse )
{     sparsity_pattern identity;
identity.resize(n, n, n);
for(size_t k = 0; k < n; k++)
identity.set(k, k, k);
fun.for_jac_sparsity(
identity, transpose, dependency, internal_bool, sparsity
);
fun.rev_hes_sparsity(
select_range, transpose, internal_bool, sparsity
);
}
else
{     b_vector select_domain(n);
for(size_t j = 0; j < n; j++)
select_domain[j] = true;
fun.for_hes_sparsity(
select_domain, select_range, internal_bool, sparsity
);
}
return;
}
// fun correspnds to g(x)
CPPAD_ASSERT_UNKNOWN( m == n );
//
// sparsity pattern for identity matrix
sparsity_pattern eye;
eye.resize(n, n, n);
for(size_t k = 0; k < n; k++)
eye.set(k, k, k);
//
if( reverse )
{     fun.rev_jac_sparsity(
eye, transpose, dependency, internal_bool, sparsity
);
}
else
{     fun.for_jac_sparsity(
eye, transpose, dependency, internal_bool, sparsity
);
}
return;
}
// ------------------------------------------------------------------------
size_t calc_hessian(
d_vector&               hessian  ,
const d_vector&         x        ,
sparse_matrix&          subset   ,
const sparsity_pattern& sparsity ,
{     size_t n_sweep;
//
if( ! global_option["hes2jac"] )
{     // fun corresponds to f(x)
//
// coloring method
std::string coloring = "cppad.symmetric";
if( global_option["colpack"] )
coloring = "colpack.symmetric";
# endif
// only one function component
d_vector w(1);
w[0] = 1.0;
//
// compute hessian
n_sweep = fun.sparse_hes(
x, w, subset, sparsity, coloring, hes_work
);
}
else
{     // fun corresponds to g(x)
//
if( global_option["subgraph"] )
{     fun.subgraph_jac_rev(x, subset);
n_sweep = 0;
}
else
{
//
// coloring method
std::string coloring = "cppad";
if( global_option["colpack"] )
coloring = "colpack";
# endif
size_t group_max = 1;
n_sweep = fun.sparse_jac_for(
group_max, x, subset, sparsity, coloring, jac_work
);
}
}
// return result
const d_vector& val( subset.val() );
size_t nnz = subset.nnz();
for(size_t k = 0; k < nnz; k++)
hessian[k] = val[k];
//
return n_sweep;
}
}

size_t                           size     ,
size_t                           repeat   ,
const CppAD::vector<size_t>&     row      ,
const CppAD::vector<size_t>&     col      ,
size_t&                          n_sweep  )
{
// --------------------------------------------------------------------
// check global options
const char* valid[] = {
"memory", "onetape", "optimize", "hes2jac", "subgraph",
"boolsparsity", "revsparsity", "subsparsity", "colpack"
# else
"boolsparsity", "revsparsity"
# endif
};
size_t n_valid = sizeof(valid) / sizeof(valid[0]);
typedef std::map<std::string, bool>::iterator iterator;
//
for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
{     if( itr->second )
{     bool ok = false;
for(size_t i = 0; i < n_valid; i++)
ok |= itr->first == valid[i];
if( ! ok )
return false;
}
}
if( global_option["subsparsity"] )
{     if( global_option["boolsparsity"] || global_option["revsparsity"] )
return false;
if( ! global_option["hes2jac"] )
return false;
}
if( global_option["subgraph"] )
{     if( ! global_option["hes2jac"] )
return false;
}
// -----------------------------------------------------------------------
// setup
size_t n = size;          // number of independent variables
CppAD::ADFun<double> fun; // AD function object used to calculate Hessian
//
// declare sparsity pattern
sparsity_pattern sparsity;
//
// declare subset where Hessian is evaluated
sparsity_pattern subset_pattern;
size_t nr  = n;
size_t nc  = n;
size_t nnz = row.size();
subset_pattern.resize(nr, nc, nnz);
for(size_t k = 0; k < nnz; k++)
subset_pattern.set(k, row[k], col[k]);
sparse_matrix subset( subset_pattern );
//
// structures that holds some of the work done by sparse_jac, sparse_hes

// -----------------------------------------------------------------------
if( ! global_option["onetape"] ) while(repeat--)
{     // choose a value for x
//
// create f(x)
create_fun(x, row, col, fun);
//
// calculate the sparsity pattern for Hessian of f(x)
calc_sparsity(sparsity, fun);
//
// calculate the Hessian at this x
jac_work.clear(); // wihtout work from previous calculation
hes_work.clear();
n_sweep = calc_hessian(
hessian, x, subset, sparsity, jac_work, hes_work, fun
);
}
else
{     // choose a value for x
//
// create f(x)
create_fun(x, row, col, fun);
//
// calculate the sparsity pattern for Hessian of f(x)
calc_sparsity(sparsity, fun);
//
while(repeat--)
{     // choose a value for x
}