The gradient flow scheme has emerged as a prominent nonperturbative renormalization scheme on the lattice, where flow time is introduced to define the renormalization scale. In this study, we perturbatively compute the gradient flow coupling for the SU(N) Yang–Mills theory in the large-N limit in terms of the lattice bare coupling up to three-loop order. This is achieved by combining the twisted Eguchi–Kawai model with the numerical stochastic perturbation theory. We analyze the flow time dependence of the perturbative coefficients to determine the perturbative beta function coefficients, successfully computing the one-loop coefficient in the large-N limit using three matrix sizes N=289,441,529. However, the higher-order coefficients are affected by large statistical errors. We also explore the potential for reducing these statistical errors through variance reduction combined with the large-N factorization property of the SU(N) Yang–Mills theory, and estimate the required number of samples for the precise determination of the higher-order coefficients.