Using finite size scaling techniques and a renormalization scheme based on the Gradient Flow, we determine non-perturbatively the $beta$-function of the $SU(3)$ Yang-Mills theory for a range of renormalized couplings $bar g^2sim 1-12$. We perform a detailed study of the matching with the asymptotic NNLO perturbative behavior at high-energy, with our non-perturbative data showing a significant deviation from the perturbative prediction down to $bar{g}^2sim1$. We conclude that schemes based on the Gradient Flow are not competitive to match with the asymptotic perturbative behavior, even when the NNLO expansion of the $beta$-function is known. On the other hand, we show that matching non-perturbatively the Gradient Flow to the Schrodinger Functional scheme allows us to make safe contact with perturbation theory with full control on truncation errors. This strategy allows us to obtain a precise determination of the $Lambda$-parameter of the $SU(3)$ Yang-Mills theory in units of a reference hadronic scale ($sqrt{8t_0},Lambda_{overline{rm MS}} = 0.6227(98)$), showing that a precision on the QCD coupling below 0.5% per-cent can be achieved using these techniques.