Prev Next Index-> contents reference index search external Up-> CppAD ADFun sparse_derivative sparse_hessian sparse_sub_hes.cpp ADFun-> record_adfun drivers Forward Reverse sparsity_pattern sparse_derivative optimize abs_normal FunCheck check_for_nan sparse_derivative-> sparse_jac sparse_jacobian sparse_hes sparse_hessian subgraph_jac_rev sparse_hessian-> sparse_hessian.cpp sub_sparse_hes.cpp sparse_sub_hes.cpp sparse_sub_hes.cpp Headings-> Purpose See Also

$\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}} }$
Subset of a Sparse Hessian: Example and Test

Purpose
This example uses a column subset of the sparsity pattern to compute a subset of the Hessian.

sub_sparse_hes.cpp
bool sparse_sub_hes(void)
{     bool ok = true;
//
// domain space vector
size_t n = 4;
for(size_t j = 0; j < n; j++)
ax[j] = double(j);

// declare independent variables and start recording

// range space vector
size_t m = 1;
ay[0] = 0.0;
for(size_t j = 0; j < n; j++)
ay[0] += double(j+1) * ax[0] * ax[j];

// create f: x -> y and stop tape recording

// sparsity pattern for the identity matrix
size_t nr     = n;
size_t nc     = n;
size_t nnz_in = n;
sparsity pattern_in(nr, nc, nnz_in);
for(size_t k = 0; k < nnz_in; k++)
{     size_t r = k;
size_t c = k;
pattern_in.set(k, r, c);
}
// compute sparsity pattern for J(x) = f'(x)
bool transpose       = false;
bool dependency      = false;
bool internal_bool   = false;
sparsity pattern_out;
f.for_jac_sparsity(
pattern_in, transpose, dependency, internal_bool, pattern_out
);
//
// compute sparsity pattern for H(x) = f''(x)
select_range[0]      = true;
f.rev_hes_sparsity(
select_range, transpose, internal_bool, pattern_out
);
size_t nnz = pattern_out.nnz();
ok        &= nnz == 7;
ok        &= pattern_out.nr() == n;
ok        &= pattern_out.nc() == n;
{     // check results
const SizeVector& row( pattern_out.row() );
const SizeVector& col( pattern_out.col() );
SizeVector row_major = pattern_out.row_major();
//
ok &= row[ row_major[0] ] ==  0  && col[ row_major[0] ] ==  0;
ok &= row[ row_major[1] ] ==  0  && col[ row_major[1] ] ==  1;
ok &= row[ row_major[2] ] ==  0  && col[ row_major[2] ] ==  2;
ok &= row[ row_major[3] ] ==  0  && col[ row_major[3] ] ==  3;
//
ok &= row[ row_major[4] ] ==  1  && col[ row_major[4] ] ==  0;
ok &= row[ row_major[5] ] ==  2  && col[ row_major[5] ] ==  0;
ok &= row[ row_major[6] ] ==  3  && col[ row_major[6] ] ==  0;
}
//
// Only interested in cross-terms. Since we are not computing rwo 0,
// we do not need sparsity entries in row 0.
for(size_t k = 0; k < 3; k++)
subset_pattern.set(k, k+1, 0);
//
// argument and weight values for computation
for(size_t j = 0; j < n; j++)
x[j] = double(n) / double(j+1);
w[0] = 1.0;
//
size_t n_sweep = f.sparse_hes(
x, w, subset, subset_pattern, coloring, work
);
ok &= n_sweep == 1;
for(size_t k = 0; k < 3; k++)
{     size_t i = k + 1;
ok &= subset.val()[k] == double(i + 1);
}
//
// convert subset from lower triangular to upper triangular
for(size_t k = 0; k < 3; k++)
subset_pattern.set(k, 0, k+1);
subset = CppAD::sparse_rcv<SizeVector, DoubleVector>( subset_pattern );
//
// This will require more work because the Hessian is computed
// column by column (not row by row).
work.clear();
n_sweep = f.sparse_hes(
x, w, subset, subset_pattern, coloring, work
);
ok &= n_sweep == 3;
//
// but it will get the right answer
for(size_t k = 0; k < 3; k++)
{     size_t i = k + 1;
ok &= subset.val()[k] == double(i + 1);
}
return ok;
}

Input File: example/sparse/sparse_sub_hes.cpp