\(\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}} }\)

valvector_llsq_obj.cpp

View page source

Using valvector to Represent Linear Least Squares Objective

# include <cppad/example/valvector/sum.hpp>
# include <cppad/example/valvector/class.hpp>
# include <cppad/cppad.hpp>
bool llsq_obj(void)
{  // ok
   bool ok = true;
   //
   // scalar_type
   typedef valvector::scalar_type scalar_type;
   //
   // eps99
   scalar_type eps99 = CppAD::numeric_limits<scalar_type>::epsilon();
   eps99            *= scalar_type(99);
   //
   // ad_valvector
   typedef CppAD::AD<valvector> ad_valvector;
   //
   // asum
   valvector_ad_sum asum;
   //
   // n_data
   // number of data points in the fit
   size_t n_data = 100;
   //
   // time, data
   // argument and data values in the function being fit
   valvector time, data;
   if( n_data == 1 )
   {  time[0] = 0.0;
      data[0] = 0.0;
   }
   else
   {  time.resize(n_data);
      data.resize(n_data);
      //
      for(size_t i = 0; i < n_data; ++i)
      {  time[i] = -1.0 + scalar_type(2 * i) / scalar_type(n_data - 1);
         if( time[i] < 0.0 )
            data[i] = -1.0;
         else if( 0.0 < time[i] )
            data[i] = +1.0;
         else
            data[i] = 0.0;
      }
   }
   // data
   if( n_data % 2 == 1 )
   {  // in case time[n_data/2] was not exactly zero due to roundoff
      data[n_data / 2] = 0.0;
   }
   //
   // nx
   // number of coefficients in the function being fit
   size_t nx = 3;
   //
   // ax
   // coefficients in the function being fit
   CPPAD_TESTVECTOR( ad_valvector ) ax(nx);
   for(size_t j = 0; j < nx; ++j)
      ax[j] = valvector(0.0);
   CppAD::Independent(ax);
   //
   // amodel
   valvector time_j(1.0);
   ad_valvector amodel(0.0);
   for(size_t  j = 0; j < nx; ++j)
   {  amodel += time_j * ax[j];
      time_j *= time;
   }
   //
   // aobj
   ad_valvector ares = data - amodel;
   ad_valvector asq  = ares * ares;
   ad_valvector aobj;
   asum(asq, aobj);
   //
   // ay
   CPPAD_TESTVECTOR( ad_valvector ) ay(1);
   ay[0] = aobj;
   CppAD::ADFun<valvector> f(ax, ay);
   //
   // x
   CPPAD_TESTVECTOR( valvector ) x(nx);
   for(size_t j = 0; j < nx; ++j)
      x[j][0] = 1.0;
   //
   // y
   CPPAD_TESTVECTOR( valvector ) y(1);
   y = f.Forward(0, x);
   //
   // res
   CPPAD_TESTVECTOR( scalar_type ) res(n_data);
   for(size_t i = 0; i < n_data; ++i)
   {  scalar_type tj = 1.0;
      scalar_type model = 0.0;
      for(size_t j = 0; j < nx; ++j)
      {  model += tj * x[j][0];
         tj    *= time[i];
      }
      res[i] = data[i] - model;
   }
   //
   // ok
   scalar_type check = 0.0;
   for(size_t i = 0; i < n_data; ++i)
      check += res[i] * res[i];
   ok   &= CppAD::NearEqual(y[0][0], check, eps99, eps99);
   //
   // dw
   CPPAD_TESTVECTOR( valvector ) w(1), dw(nx);
   w[0][0] = 1.0;
   dw = f.Reverse(1, w);
   //
   //
   // dres_sq
   CPPAD_TESTVECTOR( scalar_type ) dres_sq(nx);
   for(size_t j = 0; j < nx; ++j)
      dres_sq[j] = 0.0;
   for(size_t i = 0; i < n_data; ++i)
   {  scalar_type tj = 1.0;
      for(size_t j = 0; j < nx; ++j)
      {  dres_sq[j] += - 2.0 * res[i] * tj;
         tj         *= time[i];
      }
   }
   //
   // ok
   // reverse mode does not know that x[j] had size one and represents
   // a valvector of size n_data and with x[j][k] constant w.r.t k.
   for(size_t j = 0; j < nx; ++j)
   {  /*
      std::cout << "j = " << j;
      std::cout << ", dw[j] = " << dw[j];
      std::cout << ", dres_sq[j] = " << dres_sq[j] << "\n";
      */
      ok   &= dw[j].size() == n_data;
      ok   &= CppAD::NearEqual(dw[j].sum(), dres_sq[j], eps99, eps99);
   }
   return ok;
}