A method is given for calculating approximations to natural Riemannian cubic splines in symmetric spaces with computational efiort comparable to what is needed for the classical case of a natural cubic spline in Euclidean space. Interpolation of n+1 points in the unit sphere Sm requires the solution of a sparse linear system of 4mn linear equations. For n+1 points in bi-invariant SO(p) we have a sparse linear system of 2np(p-1) equations. Examples are given for the Euclidean sphere S2 and for bi-invariant SO(3) showing significant improvements over standard chart-based interpolants. © 2014 International Press.