We prove optimal bounds for the discretization error of geodesic finite elements for variational partial differential equations for functions that map into a nonlinear space. For this, we first generalize the well-known Céa lemma to nonlinear function spaces. In a second step, we prove optimal interpolation error estimates for pointwise interpolation by geodesic finite elements of arbitrary order. These two results are both of independent interest. Together they yield optimal a priori error estimates for a large class of manifold-valued variational problems. We measure the discretization error both intrinsically using an 𝐻1-type Finsler norm and with the 𝐻1-norm using embeddings of the codomain in a linear space. To measure the regularity of the solution, we propose a nonstandard smoothness descriptor for manifold-valued functions, which bounds additional terms not captured by Sobolev norms. As an application, we obtain optimal a priori error estimates for discretizations of smooth harmonic maps using geodesic finite elements, yielding the first high-order scheme for this problem.