  We study Hessian estimators for functions defined over an $n$-dimensional
complete analytic Riemannian manifold. We introduce new stochastic zeroth-order
Hessian estimators using $O (1)$ function evaluations. We show that, for an
analytic real-valued function $f$, our estimator achieves a bias bound of order
$ O \left( \gamma \delta^2 \right) $, where $ \gamma $ depends on both the
Levi-Civita connection and function $f$, and $\delta$ is the finite difference
step size. To the best of our knowledge, our results provide the first bias
bound for Hessian estimators that explicitly depends on the geometry of the
underlying Riemannian manifold. We also study downstream computations based on
our Hessian estimators. The supremacy of our method is evidenced by empirical
evaluations.
