Let $(M,g)$ be a Riemannian manifold and let us denote the Hessian of a function $u$ by $\mathrm{Hess}(u)$. Denote the gradient of a function $f$ by $\nabla f$. I want to prove that the formula
$2\mathrm{Hess}(u)(\nabla f, \nabla h)= g(\nabla f, \nabla g(\nabla u, \nabla h)) + g(\nabla h, \nabla g(\nabla u, \nabla f)) - g(\nabla u, \nabla g(\nabla f, \nabla h)) $
is valid for all functions $f$ and $h$ on $M$.
I am completely stuck. I have tried to compute the gradient of $g(\nabla f, \nabla h)$ for example trying to use the definition of the gradient but nothing seems to work. I would preffer not to put this in coordinates since it seems that it will be a pain, but if it is the only way then, well...
Honestly, I think index notation is less of a pain here: when we apply the product rule/metric compatiblity, it makes it easy to keep track of which slot of the Hessian we're applying the metric to. Since $f,u,h$ are all scalars, we can write all derivatives as postfix indices without getting confused, and raise and lower indices with the metric as usual. In this notation the RHS is
$$ f_i \nabla^i (u^j h_j) + h_i \nabla^i (u^j f_j) - u^j \nabla_j (f_i h^i).$$
Expanding the derivatives with the product rule gives
$$ f_i u^{ji} h_j + f^i u^j h_{ij} + h_i u_j f^{ij} + h_i f_j u^{ij} - u^j h^i f_{ij} - u_j f_i h^{ij}.$$
Cancelling the 3rd and 5th terms and likewise the 2nd and last we are left with
$$ 2 f_i h_j u^{ij} = 2 {\rm Hess}(u)(\nabla f, \nabla h)$$ as desired.
Note that we treated $f_i, u_i, h_i$ as vectors throughout - this reflects the fact that this a special case of the identity
$$ 2 (\nabla \theta) (X,Y) = g(X, \nabla (\theta(Y))) + g(Y, \nabla (\theta (X)) - \theta (\nabla (g(X,Y)))$$
with $ \theta = du, X = \nabla f, Y = \nabla g$.