Hessian matrix on Riemannian manifolds $(M, g)$ is a Riemannian manifold, and $u$ is a smooth function on $M$, one says that under a normal coordinate system, $\operatorname{Hess}(u)_{ij} = (u_{lk} - u_h \Gamma^h _{lk})B^{li}B^{kj}$, where 
$(B^{ij})$ is the square  root of matrix of $(g^{ij})$, and $(g^{ij})$ is the inverse matrix of $(g_{ij})$. How to get this conclusion?
 A: Given a smooth function $f:M\to \mathbb R$, we define on each $x\in M$ an element $df_x\in T_x^*M$ by 
$$df_x(X) = \tilde X f(x)$$
where $\tilde X$ is a vector field extended locally on some neigborhood $U$ of $x$. In local coordinate, $X = X^i \frac{\partial }{\partial x_i}$ and 
$$df_x(X)= X^i(x) \frac{\partial f}{\partial x_i}(x)$$
has nothing to do with the extension. Now we will think of $df$ as a one form on $M$, and the Hessian is defined as
$$Hf_x := \nabla df \in T_x^*M \times T_x^*M\ , \ \ \nabla df_x(X, Y) := \big( \nabla _X df\big) (Y)$$
and 
$$\big( \nabla_X df)\big) (Y):= \tilde X\big(df(\tilde Y)\big) - df\big(\nabla_X \tilde Y\big)$$ 
for some extension $\tilde X$ of $X$ and $\tilde Y$ of $Y$. In any local coordinate, 
$$Hf_{ij} = Hf\big( \frac{\partial }{\partial x_i}, \frac{\partial }{\partial x_j}\big) = 
\frac{\partial ^2 f}{\partial x_i \partial x_j} - df\big(\nabla _{\frac{\partial }{\partial x_i}} \frac{\partial }{\partial x_j}\big) = f_{ij} - \Gamma_{ij}^k f_k\ .$$
Note that your formula is only true at the center of a NORMAL coordinate as $\Gamma_{ij}^k=0$ and $g_{ij} = g^{ij} = \delta_{ij}$ at the center and thus $Hf_{ij} = f_{ij}$.
