Definition of gradient of a function $f$ in Riemannian manifold

I like to think of this in terms of matrix algebra. Let's let our vector fields and one-forms be column vectors, where the one-forms act by transpose-multiplication.

In a coordinate chart, an arbitrary vector field $V$ pulls back to a vector field on our coordinate patch in $\mathbb{R}^n$, the metric tensor pulls back to a matrix field $g$, with inverse matrix field $g^{-1}$, and the differential of $f$ is a one-form $df$. The computation defining the gradient $\nabla f$ is: $$ \langle \nabla f, V\rangle = df(V) $$ which becomes in coordinates: $$ (\nabla f)^Tg V = df^T V $$ As this is true for any $V$, we have the matrix identity $$ \nabla f^T g = df^T $$ which gives $$ g^T\nabla f = df $$ and since $g$ is a symmetric matrix, inverting we have $$ \nabla f = g^{-1}df $$

If you write it out with sums and coordinate vector fields, you will be performing this computation at the level of matrix entries, but I find this approach much cleaner.