Geometry optimization is a key aspect in any computational study of molecules. Moreover, this step may be rate-limiting in automatic potential energy surface exploration frameworks suitable for use on exascale computational resources, where potentially thousands of structures need to be optimized. Our new strategy based on geodesics dramatically speeds up optimization of molecules, thereby increasing the throughput in such applications. Optimization of molecular geometries is best done in a basis of internal coordinates: bonds lengths, bending angles, and dihedral angles. However, these coordinates are coupled, creating complicated constraints and consequently loss of efficiency for widely used optimization algorithms. Every optimization procedure works in steps that displace the molecular geometry according to the optimization goal: for minimization the step generally follows the gradients, whereas for first-order saddle point searches the step is most often based on information from the second derivate (Hessian) matrix. However, both the gradients and the Hessian (or its approximation) is only representative of the current location, and the intended step may not obey the coupling between the coordinates. This usually results in suboptimal (too short) step length or in more severe cases can lead to oscillatory behavior, eventually leading to a degraded performance. Our new optimization strategy interprets the space of all possible molecular geometries as a manifold, which is essentially the set of physically possible configurations in the space spanned by the redundant internal coordinates. The optimization steps are taken along the manifold, where the starting point is the current geometry, the end point is dictated by the optimization goal, and the path is a geodesic, i.e., the shortest path on the manifold between the starting and end points. Such a stepping strategy inherently satisfies the necessary constraints, and hence allows both better and larger steps during the optimization. Compared to the traditional approach, this substantially reduces the number of steps required to reach convergence, which we demonstrate on a molecular geometry optimization benchmark. Our method is implemented in our open-source software package, Sella. Our method can be implemented in any geometry optimization code that uses a basis of redundant internal coordinates.