Prev

Next

Index-> contents reference index search external

Up-> CppAD ADFun sparse_derivative sparse_hessian sparse_hessian.cpp

sparse_hessian.cpp

Headings

# 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;
}