SummaryThis paper addresses the problem of computing the Riemannian center of mass of a collection of symmetric positive definite matrices. We show in detail that the condition number of the Riemannian Hessian of the underlying optimization problem is never very ill conditioned in practice, which explains why the Riemannian steepest descent approach has been observed to perform well. We also show theoretically and empirically that this property is not shared by the Euclidean Hessian. We then present a limited‐memory Riemannian BFGS method to handle this computational task. We also provide methods to produce efficient numerical representations of geometric objects that are required for Riemannian optimization methods on the manifold of symmetric positive definite matrices. Through empirical results and a computational complexity analysis, we demonstrate the robust behavior of the limited‐memory Riemannian BFGS method and the efficiency of our implementation when compared to state‐of‐the‐art algorithms.