Abstract: A new approach for the effective computation of geodesic re- gression curves in shape spaces is presented. Here, one asks for a geodesic curve on the shape manifold that minimizes a sum of dissimilarity mea- sures between given two- or three-dimensional input shapes and corre- sponding shapes along the regression curve. The proposed method is based on a variational time discretization of geodesics. Curves in shape space are represented as deformations of suitable reference shapes, which renders the computation of a discrete geodesic as a PDE constrained optimization for a family of deformations. The PDE constraint is de- duced from the discretization of the covariant derivative of the velocity in the tangential direction along a geodesic. Finite elements are used for the spatial discretization, and a hierarchical minimization strategy together with a Lagrangian multiplier type gradient descent scheme is implemented. The method is applied to the analysis of root growth in botany and the morphological changes of brain structures due to aging.