Perturbative calculations of gradient flow observables are technically challenging. Current results are limited to a few quantities and, in general, to low perturbative orders. Numerical stochastic perturbation theory is a potentially powerful tool that may be applied in this context. Precise results using these techniques, however, require control over both statistical and systematic uncertainties. In this contribution, we discuss some recent algorithmic developments that lead to a substantial reduction of the cost of the computations. The matching of the ${\overline{\rm MS}}$ coupling with the gradient flow coupling in a finite box with Schr\"odinger functional boundary conditions is considered for illustration.