In this paper, a method to compute the best geodesic approximation of a set of points that belong to a Riemannian manifold is proposed. This method is based on a gradient descent technique on the tangent bundle of the manifold. An expression for the gradient is derived using the theory of Jacobi fields and an efficient numerical technique is proposed to compute these Jacobi fields. The presented approach is valid on any locally symmetric space, and the sphere S², the set of symmetric positive definite matrices P_n⁺, the special orthogonal group SO(3) and the Grassmann manifold Grass(n; p) are considered.