Finite-temperature quantum field theories are formulated in terms of Greens functions and self-energies on the Matsubara axis. In multi-orbital systems, these quantities are related to positive semidefinite matrix-valued functions of the Caratheodory and Schur class. Analysis, interpretation and evaluation of derived quantities such as real-frequency response functions requires analytic continuation of the off-diagonal elements to the real axis. We derive the criteria under which such functions exist for given Matsubara data and present an interpolation algorithm that intrinsically respects their mathematical properties. For small systems with precise Matsubara data, we find that the continuation exactly recovers all off-diagonal and diagonal elements. In real-materials systems, we show that the precision of the continuation is sufficient for the analytic continuation to commute with the Dyson equation, and we show that the commonly used truncation of off-diagonal self-energy elements leads to considerable approximation artifacts. Our method paves the way for the systematic evaluation of Matsubara data with equations of many-body theory on the real-frequency axis.