Prev Next Index-> contents reference index search external Up-> CppAD ADFun sparse_derivative sparse_hessian sparse_hessian.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_hessian.cpp Headings

$\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}} }$
Sparse Hessian: Example and Test
# include <cppad/cppad.hpp> bool sparse_hessian(void) { bool ok = true; using CppAD::AD; using CppAD::NearEqual; size_t i, j, k, ell; typedef CPPAD_TESTVECTOR(AD<double>) a_vector; typedef CPPAD_TESTVECTOR(double) d_vector; typedef CPPAD_TESTVECTOR(size_t) i_vector; typedef CPPAD_TESTVECTOR(bool) b_vector; typedef CPPAD_TESTVECTOR(std::set<size_t>) s_vector; double eps = 10. * CppAD::numeric_limits<double>::epsilon(); // domain space vector size_t n = 12; // must be greater than or equal 3; see n_sweep below 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); // range space vector size_t m = 1; a_vector a_y(m); a_y[0] = a_x[0]*a_x[1]; for(j = 0; j < n; j++) a_y[0] += a_x[j] * a_x[j] * a_x[j]; // create f: x -> y and stop tape recording // (without executing zero order forward calculation) CppAD::ADFun<double> f; f.Dependent(a_x, a_y); // new value for the independent variable vector, and weighting vector d_vector w(m), x(n); for(j = 0; j < n; j++) x[j] = double(j); w[0] = 1.0; // vector used to check the value of the hessian d_vector check(n * n); for(ell = 0; ell < n * n; ell++) check[ell] = 0.0; ell = 0 * n + 1; check[ell] = 1.0; ell = 1 * n + 0; check[ell] = 1.0 ; for(j = 0; j < n; j++) { ell = j * n + j; check[ell] = 6.0 * x[j]; } // ------------------------------------------------------------------- // second derivative of y[0] w.r.t x d_vector hes(n * n); hes = f.SparseHessian(x, w); for(ell = 0; ell < n * n; ell++) ok &= NearEqual(w[0] * check[ell], hes[ell], eps, eps ); // -------------------------------------------------------------------- // example using vectors of bools to compute sparsity pattern for Hessian b_vector r_bool(n * n); for(i = 0; i < n; i++) { for(j = 0; j < n; j++) r_bool[i * n + j] = false; r_bool[i * n + i] = true; } f.ForSparseJac(n, r_bool); // b_vector s_bool(m); for(i = 0; i < m; i++) s_bool[i] = w[i] != 0; b_vector p_bool = f.RevSparseHes(n, s_bool); hes = f.SparseHessian(x, w, p_bool); for(ell = 0; ell < n * n; ell++) ok &= NearEqual(w[0] * check[ell], hes[ell], eps, eps ); // -------------------------------------------------------------------- // example using vectors of sets to compute sparsity pattern for Hessian s_vector r_set(n); for(i = 0; i < n; i++) r_set[i].insert(i); f.ForSparseJac(n, r_set); // s_vector s_set(m); for(i = 0; i < m; i++) if( w[i] != 0. ) s_set[0].insert(i); s_vector p_set = f.RevSparseHes(n, s_set); // example passing sparsity pattern to SparseHessian hes = f.SparseHessian(x, w, p_set); for(ell = 0; ell < n * n; ell++) ok &= NearEqual(w[0] * check[ell], hes[ell], eps, eps ); // -------------------------------------------------------------------- // use row and column indices to specify upper triangle of // non-zero elements of Hessian size_t K = n + 1; i_vector row(K), col(K); hes.resize(K); k = 0; for(j = 0; j < n; j++) { // diagonal of Hessian row[k] = j; col[k] = j; k++; } // only off diagonal non-zero elemenet in upper triangle row[k] = 0; col[k] = 1; k++; ok &= k == K; CppAD::sparse_hessian_work work; // can use p_set or p_bool. size_t n_sweep = f.SparseHessian(x, w, p_set, row, col, hes, work); for(k = 0; k < K; k++) { ell = row[k] * n + col[k]; ok &= NearEqual(w[0] * check[ell], hes[k], eps, eps ); } ok &= n_sweep == 2; // now recompute at a different x and w (using work from previous call w[0] = 2.0; x[1] = 0.5; ell = 1 * n + 1; check[ell] = 6.0 * x[1]; s_vector not_used; n_sweep = f.SparseHessian(x, w, not_used, row, col, hes, work); for(k = 0; k < K; k++) { ell = row[k] * n + col[k]; ok &= NearEqual(w[0] * check[ell], hes[k], eps, eps ); } ok &= n_sweep == 2; return ok; } 
Input File: example/sparse/sparse_hessian.cpp