calculus and analysis - Compute Hessian of function symbolically

This question could equally apply to the computation of other symbolic transformations of a function, but I use the Hessian as an example here.

Consider a two-variable function fun which could be defined either using down-values or as a Function. It may or may not contain symbolic parameters.

Compute the Hessian of fun and return it as a Function. The symbolic computation must already be done and the result embedded in the result function. This function will be used for multiple numerical evaluations later.

---

One possible way to do this is:

hessian[fun_] := 
  Block[{hess, x, y}, 
    Function[{x, y}, hess] /. hess -> D[fun[x, y], {{x, y}, 2}]
  ]

What's wrong with this? Consider cases when fun may contain parameters (which in my case will usually be global variables with their value already set). What if these conflict with the Block variables?

x = 5;
hessian[Function[{a, b}, x a^3 + b^3]]
(* Function[{x, y}, {{12 x^2, 0}, {0, 6 y}}] *)


p = 5;
hessian[Function[{a, b}, p a^3 + b^3]]
(* Function[{x, y}, {{30 x, 0}, {0, 6 y}}] *)

One workaround is to make these live in a context that is not user-accessible, e.g.

hessian[fun_] := 
 Block[{hessian`hess, hessian`x, hessian`y}, 
  Function[{hessian`x, hessian`y}, hessian`hess] /. 
   hessian`hess -> 

    D[fun[hessian`x, hessian`y], {{hessian`x, hessian`y}, 2}]]

But this is getting ugly quickly. What is a better way?

Answer

I think this should be fine, the full form doesn't look neat with x$123 etc but that doesn't matter:

hessian[fun_] := Module[{x, y}, 
    Function @@ {{x, y}, D[fun[x, y], {{x, y}, 2}]}
]

this will work too:

 ... Function[{x, y}, Evaluate @ D[fun[x, y], {{x, y}, 2}]] ...

and will work even nicer because the outer Module will rename Function's variables to x$ without numbers.