[ADMB Users] how to calculate the hessian matrix by autodif library

[dave fournier]

sorry the entire example did not paste here it is again


#include <admodel.h>


dmatrix fmm::hessian(void)
{
   dmatrix tmp(1,n,1,n);
   for (int i=1;i<=n;i++)
     for (int j=1;j<=n;j++)
       tmp(i,j)=h(i,j);
   return tmp;
}

double func2(const dvar_vector & x)
{
   dvariable f;
   f=norm2(x-1.0);
   return value(f);
}

int main()
{
     int nvar=5;
     independent_variables x(1,nvar);
     x.initialize();
     double f;
     dvector g(1,nvar);
     dmatrix dHessian(1,nvar,1,nvar);
     fmm fmc(nvar); // creates the function minimizing control structure
     fmc.crit=0.000001;
   /////////////////////////////////////////////////////
   gradient_structure gs;
    while(fmc.ireturn>=0)//
   {//
    fmc.fmin(f,x,g);

   if(fmc.ireturn>0)
    {
            f=func2(x);
      gradcalc(nvar,g);
        }
   }
   dHessian=fmc.hessian();
   cout << "The minimum value = " << "at x =\n"
<< x <<"\n"<<"G"<<g<<f<< "\n";
   cout << "fmm estimate of Hessian" << endl;
   cout << setw(10) << dHessian << endl;
   dvector g1(1,nvar);
   dvector xsave(1,nvar);
   xsave=x;
   for (int i=1;i<=nvar;i++)
   {
     x=xsave;
     double delta=1.e-5;
     x(i)+=delta;
     f=func2(x);
     gradcalc(nvar,g);
     x=xsave;
     x(i)-=delta;
     f=func2(x);
     gradcalc(nvar,g1);
     dHessian(i)=(g-g1)/(2.0*delta);
   }
   cout << "finite difference estimate of Hessian" << endl;
   cout << setw(10) << dHessian << endl;

   return 0;
}