Prev Next repeat_newton_hessian

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@
Repeated Computation of Control Problem Hessian Using Newton Method

Syntax
repeat_kedem_hessian(repeat, J, full_step, x, hess)
 repeat_partial_hessian(repeat, size)
 repeat_full_hessian(repeat, size)
repeat
is the number of times to choose a new value for the controls.

J
is the number of time points in the discrete version of the control problem.

full_step
If this is true, the full Newton step is used to compute the Hessian. Otherwise, the partial Newton step mode is used.

x
If x has size zero, it is not used. Otherwise, it must have size 2*J . In this case it specifies the value of the control vector for the first Hessian computation.

hess
This vector must have size 4*J*J . The input value of its elements does not matter. Upon return, it contains the value last Hessian computed by repeat_newton_hessian.

size
This is the number time points in the problem; i.e., the value of J for this test. 
void repeat_newton_hessian(
     size_t                repeat     ,
     size_t                J          ,
     bool                  full_step     ,
     const VECTOR(double)& x          ,
     VECTOR(double)&       hess       )
{    assert( x.size() == 0 || size_t( x.size() ) == 2 * J );
     assert( size_t( hess.size() ) == 2 * J * 2 * J );
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
     typedef CppAD::AD<double> adouble;
     // ----------------------------------------------------------------------
     // control problem parameters
     double T;
     VECTOR(double) p(4), q(4);
     set_T_p_and_q(T, p, q);
     double delta_t = T / double(J - 1);
     // ----------------------------------------------------------------------
     // L_fun, F_fun
     CppAD::ADFun<double> F_fun, L_fun;
     control::rec_objective(F_fun, J,  delta_t, q);
     control::rec_constraint(L_fun, J, delta_t, p);
     size_t m = L_fun.Range();
     size_t n = L_fun.Domain() - m;
     assert( F_fun.Range() == 1 );
     assert( F_fun.Domain() == n + m );
     // -----------------------------------------------------------------------
     // aL_fun
     CppAD::ADFun<adouble> aL_fun;
     adouble adelta_t = delta_t;
     VECTOR(adouble) ap(4);
     for(size_t i = 0; i < 4; i++)
          ap[i] = p[i];
     control::rec_constraint(aL_fun, J, adelta_t, ap);
     // -----------------------------------------------------------------------
     // control_solve
     double criteria = double(n + m) * eps99;
     control::implicit_solver control_solve(L_fun, aL_fun, criteria);
     // -----------------------------------------------------------------------
     // R_newton
     size_t num_step = 2;
     implicit_newton<control::implicit_solver> R_newton(
          full_step, num_step, aL_fun, F_fun, control_solve
     );
     //
     VECTOR(double) x0(n), x1(n), w(2), dRR(2 * n);
     w[0] = 0.0;
     w[1] = 1.0;
     for(size_t i = 0; i < n; i++)
          x1[i] = 0.0;
     if( x.size() == 0 )
          CppAD::uniform_01(n, x0);
     else
          x0 = x;
     for(size_t count = 0; count < repeat; count++)
     {    // zero order forward
          R_newton.Forward(0, x0);
          for(size_t i = 0; i < n; i++)
          {    x1[i] = 1.0;
               R_newton.Forward(1, x1);
               dRR = R_newton.Reverse(2, w);
               for(size_t j = 0; j < n; j++)
                    hess[i * n + j] = dRR[2 * j + 0];
               x1[i] = 0.0;
          }
          // pick a random value for next x0
          CppAD::uniform_01(n, x0);
     }
}
void time_partial_hessian(size_t size, size_t repeat)
{    size_t J = size;
     size_t full_step = false;
     VECTOR(double) hess(2 * J * 2 * J), x(0);
     repeat_newton_hessian(repeat, J, full_step, x, hess);
}
void time_full_hessian(size_t size, size_t repeat)
{    size_t J = size;
     size_t full_step = true;
     VECTOR(double) hess(2 * J * 2 * J), x(0);
     repeat_newton_hessian(repeat, J, full_step, x, hess);
}
Input File: src/time.cpp