Efficient computation of a single high-order partial derivative

[kasselhingee]

When computing the probability density function $f(u)$ from a cumulative density function $F(u)$ a single high-order partial derivative is needed. When $u$ is a vector of length $d$, the density is:

$$f(u) = \frac{\partial^d} { \partial u_1 \partial u_2 \dots \partial u_d} F(u)$$

My colleague and I want to avoid computing all $d$-order partial derivatives because there are a lot of them and we only want one.

The best idea I have so far is:

1. apply Forward(1, e_1) following this help page where $e_1$ is the first basis vector: (1, 0, 0, .... 0).
2. optimize
3. base2ad
4. repeat for e_2, e_3, ... e_d.

I'm wondering if this is a good way to go about this, and it there are other ways to go about this.

[bradbell]

Here is one possible improvement. The idea is that each derivative only has one variable and n-1 dynamic parameters. This may greatly reduce the amount of memory for the derivative computation. On the other hand, this method has twice as many base2ad operations, so it may not be an improvement:

# include <cppad/cppad.hpp>

namespace {
   // vector, AD
   using CppAD::vector;
   using CppAD::AD;
   //
   // prod
   AD<double> prod(const vector< AD<double> >& au)
   {  // n
      size_t n = au.size();
      //
      // result
      AD<double> result = 1.0;
      for(size_t j = 0; j < n; ++j)
         result *= au[j];
      //
      return result;
   }
   //
   // combine
   template <class Vector>
   void combine(Vector&u, const Vector& x, const Vector& p, size_t x_index )
   {  assert( x.size() == 1 );
      assert( u.size() == p.size() + 1 );
      assert( x_index < u.size() );
      //
      // n
      size_t n = u.size();
      //
      // u
      for(size_t j = 0; j < x_index; ++j)
         u[j] = p[j];
      u[x_index] = x[0];
      for(size_t j = x_index + 1; j < n; ++j)
         u[j] = p[j-1];
      //
      return;
   }
   //
   // split
   template <class Vector>
   void split(const Vector&u, Vector& x, Vector& p, size_t x_index )
   {  assert( x.size() == 1 );
      assert( u.size() == p.size() + 1 );
      assert( x_index < u.size() );
      //
      // n
      size_t n = u.size();
      //
      // u
      for(size_t j = 0; j < x_index; ++j)
         p[j] = u[j];
      x[0] = u[x_index];
      for(size_t j = x_index + 1; j < n; ++j)
         p[j-1] = u[j];
      //
      return;
   }
}

int main(void)
{  // n
   size_t n = 4;
   //
   // au, ax, ap, abort_op_index, record_compare
   vector< AD<double> > au(n), ax(1), ap(n-1);
   for(size_t j = 0; j < n; ++j)
      au[j] = double(j + 1);
   split(au, ax, ap, 0);
   size_t abort_op_index = 0;
   bool   record_compare = false;
   //
   // ae
   vector< AD<double> > ae(1);
   ae[0] = 1.0;
   //
   // ay
   vector< AD<double> > ay(1);
   //
   // ax, ap
   CppAD::Independent(ax, abort_op_index, record_compare, ap);
   //
   // f = prod(x[0], p[0], p[1] ...) 
   combine(au, ax, ap, 0);
   ay[0] = prod(au);
   CppAD::ADFun<double> f(ax, ay);
   //
   // f = patial_x[0] prod(x[0], p[0], p[1] ...) 
   CppAD::ADFun< AD<double> , double>  af = f.base2ad();
   CppAD::Independent(ax, abort_op_index, record_compare, ap);
   af.new_dynamic(ap);
   af.Forward(0, ax);
   ay = af.Forward(1, ae);
   f.Dependent(ax, ay);
   //
   // ax_old, ap_old
   vector< AD<double> > ax_old(1), ap_old(n-1);
   //
   // f = patial_p[0]  prod(p[0], x[0], p[1] ...) 
   af = f.base2ad();
   CppAD::Independent(ax, abort_op_index, record_compare, ap);
   ax_old[0] = ap[0];
   ap_old[0] = ax[0];
   for(size_t j = 1; j < n-1; ++j)
      ap_old[j] = ap[j];
   af.new_dynamic(ap_old);
   ay = af.Forward(0, ax_old);
   f.Dependent(ax, ay);
   //
   // f = partial_p[0] patial_x[0] prod(p[0], x[0], p[1] ...) 
   af = f.base2ad();
   CppAD::Independent(ax, abort_op_index, record_compare, ap);
   af.new_dynamic(ap);
   af.Forward(0, ax);
   ay = af.Forward(1, ae);
   f.Dependent(ax, ay);
   //
   // result
   vector<double> result;
   vector<double> u(n), x(1), p(n-1);
   for(size_t j = 0; j < n; ++j)
      u[j] = double(j + 1);
   split(u, x, p, 1);
   f.new_dynamic(p);
   result = f.Forward(0, x);
   //
   // check
   // start j at 2 because we have taken the partial w.r.t. u[0] and u[1]
   double check = 1.0; 
   for(size_t j = 2; j < n; ++j)
      check *= u[j]; 
   //
   std::cout << "result = " << result[0] << "\n";
   std::cout << "check  = " << check << "\n";
   //
   return 0;
}
cppad>

[kasselhingee]

Thank you @bradbell ! I like your recording of a new function with independent and dynamic parameters swapped. I think it will work well as I'm not too worried about base2ad operations to obtain the final function, more about the operations required to evaluate the final function.