  In this article, we propose and analyze an effective Hessian recovery
strategy for the Lagrangian finite element of arbitrary order $k$. We prove
that the proposed Hessian recovery preserves polynomials of degree $k+1$ on
general unstructured meshes and superconverges at rate $O(h^k)$ on mildly
structured meshes. In addition, the method preserves polynomials of degree
$k+2$ on translation invariant meshes and produces a symmetric Hessian matrix
when the sampling points for recovery are selected with symmetry. Numerical
examples are presented to support our theoretical results.
