\(\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_sparse_hessian.cpp¶
View page sourceCppad Speed: Sparse Hessian¶
Specifications¶
See link_sparse_hessian .
Implementation¶
# include <cppad/cppad.hpp>
# include <cppad/speed/uniform_01.hpp>
# include <cppad/speed/sparse_hes_fun.hpp>
// Note that CppAD uses global_option["memory"] at the main program level
# include <map>
extern std::map<std::string, bool> global_option;
// see comments in main program for this external
extern size_t global_cppad_thread_alloc_inuse;
namespace {
// typedefs
using CppAD::vector;
typedef CppAD::AD<double> a1double;
typedef CppAD::AD<a1double> a2double;
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_rc<s_vector> sparsity_pattern;
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 ,
CppAD::ADFun<double>& fun )
{
// 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_conditional_skip no_compare_op no_print_for_op";
if( global_option["val_graph"] )
optimize_options += " val_graph";
//
// order of derivative in sparse_hes_fun
size_t order = 0;
//
// do not even record comparison operators
size_t abort_op_index = 0;
bool record_compare = false;
//
if( ! global_option["hes2jac"] )
{
// declare independent variables
Independent(a1x, abort_op_index, record_compare);
//
// 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, abort_op_index, record_compare);
//
// 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)
CppAD::ADFun<a1double> a1f;
a1f.Dependent(a2x, a2y);
//
// declare independent variables for g(x)
Independent(a1x, abort_op_index, record_compare);
//
// 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 ,
CppAD::ADFun<double>& fun )
{
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 ,
CppAD::sparse_jac_work& jac_work ,
CppAD::sparse_hes_work& hes_work ,
CppAD::ADFun<double>& fun )
{ size_t n_color;
//
if( ! global_option["hes2jac"] )
{ // fun corresponds to f(x)
//
// coloring method
std::string coloring = "cppad";
if( global_option["colpack"] )
coloring = "colpack";
if( global_option["symmetric"] )
coloring += ".symmetric";
else
coloring += ".general";
//
// only one function component
d_vector w(1);
w[0] = 1.0;
//
// compute hessian
n_color = 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_color = 0;
}
else
{
//
// coloring method
std::string coloring = "cppad";
# if CPPAD_HAS_COLPACK
if( global_option["colpack"] )
coloring = "colpack";
# endif
size_t group_max = 1;
n_color = 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_color;
}
}
bool link_sparse_hessian(
size_t size ,
size_t repeat ,
const CppAD::vector<size_t>& row ,
const CppAD::vector<size_t>& col ,
CppAD::vector<double>& x ,
CppAD::vector<double>& hessian ,
size_t& n_color )
{ global_cppad_thread_alloc_inuse = 0;
// --------------------------------------------------------------------
// check global options
const char* valid[] = {
"memory", "onetape", "optimize", "hes2jac", "subgraph",
"boolsparsity", "revsparsity", "symmetric", "val_graph"
# if CPPAD_HAS_COLPACK
, "colpack"
# else
, "subsparsity"
# 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;
}
# if ! CPPAD_HAS_COLPACK
if( global_option["colpack"] )
return false;
# endif
// -----------------------------------------------------------------------
// 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
CppAD::sparse_jac_work jac_work;
CppAD::sparse_hes_work hes_work;
// -----------------------------------------------------------------------
if( ! global_option["onetape"] ) while(repeat--)
{ // choose a value for x
CppAD::uniform_01(n, 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_color = calc_hessian(
hessian, x, subset, sparsity, jac_work, hes_work, fun
);
}
else
{ // choose a value for x
CppAD::uniform_01(n, 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
CppAD::uniform_01(n, x);
//
// calculate this Hessian at this x
n_color = calc_hessian(
hessian, x, subset, sparsity, jac_work, hes_work, fun
);
}
}
size_t thread = CppAD::thread_alloc::thread_num();
global_cppad_thread_alloc_inuse = CppAD::thread_alloc::inuse(thread);
return true;
}