It sounds like you want the Hessian of a Riemannian manifold.
(If you just want the Hessian in spherical coordinates in $\mathbb{R}^3$, then you should just use the chain rule on the regular Hessian).
In this case, the Hessian is written: $$ \mathcal{H}[f] = (\partial_{ij}f-\Gamma^k_{ij}\partial_kf)[dx^i\otimes dx^j] $$
In local coordinates, the unit sphere is given by: $$ \Psi(\theta,\phi)=(\cos(\theta)\sin(\phi),\sin(\theta)\sin(\phi),\cos(\phi)) $$
So the metric tensor is written: $$ g=\begin{bmatrix}\sin^2(\phi) & 0 \\ 0 & 1\end{bmatrix} $$
with Christoffel symbols given by (e.g. see here):
\begin{align}
\Gamma_{11}^1 &= \Gamma_{22}^1 = \Gamma_{12}^2 = \Gamma_{22}^2=0\\[2mm]
\Gamma_{12}^1 &= \frac{\cos(\phi)}{\sin(\phi)}\\[1mm]
\Gamma_{11}^2 &= -\sin(\phi)\cos(\phi)
\end{align}
which gives me the result:
$$
\mathcal{H}[f] = \begin{bmatrix}
\partial_{\theta\theta}f + \cos(\phi)\sin(\phi)\partial_\phi f & \;\partial_{\theta\phi}f-\dfrac{\cos(\phi)}{\sin(\phi)}\partial_\theta f\\
\partial_{\theta\phi}f-\dfrac{\cos(\phi)}{\sin(\phi)}\partial_\theta f & \;\partial_{\phi\phi}f
\end{bmatrix}
$$
Apologies if there are any computation errors.
See also: [1], [2].
One check is that the Laplace-Beltrami operator is the trace of the Hessian, which is written as
$$ \text{tr}(\mathcal{H}[f]) = g^{ij}\mathcal{H}[f]_{ij} =
\frac{1}{\sin^2(\phi)}\left[ \partial_{\theta\theta}f+\sin(\phi)\cos(\phi)\partial_\phi f \right] + \partial_{\phi\phi}f $$
From here, it is given by:
$$
\Delta_g f = (\sin(\phi))^{-1}{\partial}_{\phi}\left(
\sin(\phi)\partial_\phi f
\right) + (\sin(\phi))^{-2}\partial_{\theta\theta}f
$$
which corroborates it.
Compute the eigenvalues of the hessian.
If all the eigenvalues are nonnegative, it is positive semidefinite.
If all the eigenvalues are positive, it is positive definite.
If all the eigenvalues are nonpositive, it is negative semidefinite.
If all the eigenvalues are negative, it is negative definite.
Otherwise, it is indefinite.
Edit:
For that example, you have found $c=(0,0,0)$.
$$H(f(c))=\begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0\end{bmatrix}$$
$$H(f(c))-\lambda I= \begin{bmatrix} -\lambda & 1 & 1 \\ 1 & -\lambda & 0 \\ 1 & 0 & -\lambda\end{bmatrix}$$
\begin{align}\det(H(f(c))-\lambda I)&= \det\left(\begin{bmatrix}1 & 1 \\ -\lambda & 0 \end{bmatrix} \right) -\lambda \det \left( \begin{bmatrix}-\lambda & 1 \\ 1& -\lambda \end{bmatrix} \right)
\\&=\lambda-\lambda(\lambda^2-1)
\\&=\lambda(2-\lambda^2)\end{align}
Hence, the eigenvalues are $0$, $\sqrt{2}$ and $-\sqrt{2}$. Hence it is indefinite.
Best Answer
The statement should probably be viewed in the context of the book. When e.g. the matrix is diagonally dominated you may use the diagonal inverse to calculate the inverse using a von Neumann series. So it is 'sort' of an approximation. Better if the diagonal elements dominate with some large factor the other elements. But in general, you are right that the formula need not make sense.