We recover the Riemannian gradient of a given function defined on interior points of a Riemannian submanifold in the Euclidean space based on a sample of function evaluations at points in the submanifold. This approach is based on the estimates of the Laplace-Beltrami operator proposed in the diffusion-maps theory. The Riemannian gradient estimates do not involve differential terms. Analytical convergence results of the Riemannian gradient expansion are proved. We apply the Riemannian gradient estimate in a gradient-based algorithm providing a derivative-free optimization method. We test and validate several applications, including tomographic reconstruction from an unknown random angle distribution, and the sphere packing problem in dimens...