Matching on covariates is a well-established framework for estimating causal effects in observational studies. The principal challenge stems from the often high-dimensional structure of the problem. Many methods have been introduced to address this, with different advantages and drawbacks in computational and statistical performance as well as interpretability. This article introduces a natural optimal matching method based on multimarginal unbalanced optimal transport that possesses many useful properties in this regard. It provides interpretable weights based on the distance of matched individuals, can be efficiently implemented via the iterative proportional fitting procedure, and can match several treatment arms simultaneously. Importantly, the proposed method only selects good matches from either group, hence is competitive with the classical k-nearest neighbors approach in terms of bias and variance in finite samples. Moreover, we prove a central limit theorem for the empirical process of the potential functions of the optimal coupling in the unbalanced optimal transport problem with a fixed penalty term. This implies a parametric rate of convergence of the empirically obtained weights to the optimal weights in the population for a fixed penalty term.