title: How to compile a function that computes the Hessian?
body: <p>I am looking to see how a function that computes the Hessian of a log-likelihood can be compiled, so that it can be efficiently used with different sets of parameters. </p>
  
  <p>Here is an example. </p>
  
  <p>Suppose we have a function that computes the log-likelihood of a logit model, where y is vector and x is a matrix. beta is a vector of parameters.  </p>
  
  <pre><code> pLike[y_, x_, beta_] :=
   Module[
    {xbeta, logDen},
    xbeta = x.beta;
    logDen = Log[1.0 + Exp[xbeta]];
    Total[y*xbeta - logDen]
    ]
  </code></pre>
  
  <p>Given the following data, we can use it as follows</p>
  
  <pre><code>In[1]:= beta = {0.5, -1.0, 1.0};
  
  In[2]:= xmat = 
    Table[Flatten[{1, 
       RandomVariate[NormalDistribution[0.0, 1.0], {2}]}], {500}];
  
  In[3]:= xbeta = xmat.beta;
  
  In[4]:= prob = Exp[xbeta]/(1.0 + Exp[xbeta]);
  
  In[5]:= y = Map[RandomVariate[BernoulliDistribution[#]] &amp;, prob] ;
  
  In[6]:= Tally[y]
  
  Out[6]= {{1, 313}, {0, 187}}
  
  In[9]:= pLike[y, xmat, beta]
  
  Out[9]= -272.721
  </code></pre>
  
  <p>We can write its hessian as follows</p>
  
  <pre><code> hessian[y_, x_, z_] :=
    Module[{},
     D[pLike[y, x, z], {z, 2}]
    ]
  
  
  In[10]:= z = {z1, z2, z3}
  
  Out[10]= {z1, z2, z3}
  
  In[11]:= AbsoluteTiming[hess = hessian[y, xmat, z];]
  
  Out[11]= {0.1248040, Null}
  
  In[12]:= AbsoluteTiming[
   Table[hess /. {z1 -&gt; 0.0, z2 -&gt; -0.5, z3 -&gt; 0.8}, {100}];]
  
  Out[12]= {14.3524600, Null}
  </code></pre>
  
  <p>For efficiency reasons, I can compile the original likelihood function as follows</p>
  
  <pre><code>pLikeC = Compile[{{y, _Real, 1}, {x, _Real, 2}, {beta, _Real, 1}},
     Module[
      {xbeta, logDen},
      xbeta = x.beta;
      logDen = Log[1.0 + Exp[xbeta]];
      Total[y*xbeta - logDen]
      ],
     CompilationTarget -&gt; "C", Parallelization -&gt; True,  
     RuntimeAttributes -&gt; {Listable}
     ];
  </code></pre>
  
  <p>which yields the same answer as pLike</p>
  
  <pre><code>In[10]:= pLikeC[y, xmat, beta]
  
  Out[10]= -272.721
  </code></pre>
  
  <p>I am looking for an easy way to obtain similarly, a compiled version of the hessian function, given my interest in evaluating it many times. </p>