The non-abelian Yang-Mills self-interaction is marginal in four dimensions in the usual terminology of the renormalization group, and thus requires renormalization to all orders in an expansion (in either the coupling or a loop expansion). It is well-known in principle how to reconcile UV-renormalization with local gauge invariance, but a fully rigorous and satisfactory treatment taking also into account the problems posed by the complicated IR-behavior of the theory on non-compact manifolds has only been given relatively recently. In this talk, I describe how a combination of the Epstein-Glaser method/method of renormalization group flow equations and cohomological techniques a la Batalin-Vilkovisky can be used to establish the existence of correlation n-point functions of arbitrary gauge invariant local operators to arbitrary orders in perturbation theory. The second method of proof comes by construction with rather refined bounds for these objects capturing their short and long distance properties. Time permitting, I furthermore explain how the operator product arises in this context.