The theory of geodesic regression aims to find a geodesic curve which is an optimal fit to a given set of data. In this article we restrict ourselves to the Riemannian manifold of positive definite operators (matrices) on a Hilbert space of finite dimension. There is a unique geodesic curve connecting two positive definite operators, and it is given by the weighted geometric mean. The function that measures the squared Riemannian metric distance between an operator and a geodesic curve is not convex nor geodesically convex in the operators generating the curve. This is a marked difference to the situation in linear regression. The literature mainly tries to find numerical solutions that approximate the optimal curve in a single point. We suggest to apply a distance measure slightly coarser than the Riemannian metric. The ensuing control function faithfully identifies geodesic curves, and it coincides with the standard control function based on the Riemannian metric for commuting operators. The control function constructed in this way is geodesically convex. We are therefore able to find a global and uniquely defined optimal fit to any given set of data. The generators of the geodesic curve may also be determined as the unique solution to two operator equations.