Prev Next bender_quad.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}} }@)@ This is cppad-20221105 documentation. Here is a link to its current documentation .
BenderQuad: Example and Test
Define @(@ F : \B{R} \times \B{R} \rightarrow \B{R} @)@ by @[@ F(x, y) = \frac{1}{2} \sum_{i=1}^N [ y * \sin ( x * t_i ) - z_i ]^2 @]@ where @(@ z \in \B{R}^N @)@ is a fixed vector. It follows that @[@ \begin{array}{rcl} \partial_y F(x, y) & = & \sum_{i=1}^N [ y * \sin ( x * t_i ) - z_i ] \sin( x * t_i ) \\ \partial_y \partial_y F(x, y) & = & \sum_{i=1}^N \sin ( x t_i )^2 \end{array} @]@ Furthermore if we define @(@ Y(x) @)@ as the argmin of @(@ F(x, y) @)@ with respect to @(@ y @)@, @[@ \begin{array}{rcl} Y(x) & = & y - [ \partial_y \partial_y F(x, y) ]^{-1} \partial_y F[x, y] \\ & = & \left. \sum_{i=1}^N z_i \sin ( x t_i ) \right/ \sum_{i=1}^N z_i \sin ( x * t_i )^2 \end{array} @]@ 

# include <cppad/cppad.hpp>

namespace {
    using CppAD::AD;
    typedef CPPAD_TESTVECTOR(double)         BAvector;
    typedef CPPAD_TESTVECTOR(AD<double>)   ADvector;

    class Fun {
    private:
        BAvector t_; // measurement times
        BAvector z_; // measurement values
    public:
        // constructor
        Fun(const BAvector &t, const BAvector &z)
        : t_(t), z_(z)
        { }
        // Fun.f(x, y) = F(x, y)
        ADvector f(const ADvector &x, const ADvector &y)
        {   size_t i;
            size_t N = size_t(z_.size());

            ADvector F(1);
            F[0] = 0.;

            AD<double> residual;
            for(i = 0; i < N; i++)
            {   residual = y[0] * sin( x[0] * t_[i] ) - z_[i];
                F[0]    += .5 * residual * residual;
            }
            return F;
        }
        // Fun.h(x, y) = H(x, y) = F_y (x, y)
        ADvector h(const ADvector &x, const BAvector &y)
        {   size_t i;
            size_t N = size_t(z_.size());

            ADvector fy(1);
            fy[0] = 0.;

            AD<double> residual;
            for(i = 0; i < N; i++)
            {   residual = y[0] * sin( x[0] * t_[i] ) - z_[i];
                fy[0]   += residual * sin( x[0] * t_[i] );
            }
            return fy;
        }
        // Fun.dy(x, y, h) = - H_y (x,y)^{-1} * h
        //                 = - F_yy (x, y)^{-1} * h
        ADvector dy(
            const BAvector &x ,
            const BAvector &y ,
            const ADvector &H )
        {   size_t i;
            size_t N = size_t(z_.size());

            ADvector Dy(1);
            AD<double> fyy = 0.;

            for(i = 0; i < N; i++)
            {   fyy += sin( x[0] * t_[i] ) * sin( x[0] * t_[i] );
            }
            Dy[0] = - H[0] / fyy;

            return Dy;
        }
    };

    // Used to test calculation of Hessian of G
    AD<double> G(const ADvector& x, const BAvector& t, const BAvector& z)
    {   // compute Y(x)
        AD<double> numerator = 0.;
        AD<double> denominator = 0.;
        size_t k;
        for(k = 0; k < size_t(t.size()); k++)
        {   numerator   += sin( x[0] * t[k] ) * z[k];
            denominator += sin( x[0] * t[k] ) * sin( x[0] * t[k] );
        }
        AD<double> y = numerator / denominator;

        // V(x) = F[x, Y(x)]
        AD<double> sum = 0;
        for(k = 0; k < size_t(t.size()); k++)
        {   AD<double> residual = y * sin( x[0] * t[k] ) - z[k];
            sum += .5 * residual * residual;
        }
        return sum;
    }
}

bool BenderQuad(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;

    // temporary indices
    size_t i, j;

    // x space vector
    size_t n = 1;
    BAvector x(n);
    x[0] = 2. * 3.141592653;

    // y space vector
    size_t m = 1;
    BAvector y(m);
    y[0] = 1.;

    // t and z vectors
    size_t N = 10;
    BAvector t(N);
    BAvector z(N);
    for(i = 0; i < N; i++)
    {   t[i] = double(i) / double(N);       // time of measurement
        z[i] = y[0] * sin( x[0] * t[i] );   // data without noise
    }

    // construct the function object
    Fun fun(t, z);

    // evaluate the G(x), G'(x) and G''(x)
    BAvector g(1), gx(n), gxx(n * n);
    CppAD::BenderQuad(x, y, fun, g, gx, gxx);


    // create ADFun object Gfun corresponding to G(x)
    ADvector a_x(n), a_g(1);
    for(j = 0; j < n; j++)
        a_x[j] = x[j];
    Independent(a_x);
    a_g[0] = G(a_x, t, z);
    CppAD::ADFun<double> Gfun(a_x, a_g);

    // accuracy for checks
    double eps = 10. * CppAD::numeric_limits<double>::epsilon();

    // check Jacobian
    BAvector check_gx = Gfun.Jacobian(x);
    for(j = 0; j < n; j++)
        ok &= NearEqual(gx[j], check_gx[j], eps, eps);

    // check Hessian
    BAvector check_gxx = Gfun.Hessian(x, 0);
    for(j = 0; j < n*n; j++)
        ok &= NearEqual(gxx[j], check_gxx[j], eps, eps);

    return ok;
}
Input File: example/general/bender_quad.cpp