The gradient flow method is a renormalization scheme in which the gauge field is flowed by the diffusion equation.
The gradient flow scheme has benefits that the observables composed of flowed gauge fields
do not require further renormalization and do not depend on the regularization.
From the independence of the regularization, this scheme allows us to relate
the lattice regularization and the dimensional regularization such as the \MSbar scheme.
We compute the gradient flow coupling for the twisted Eguchi--Kawai model
using the numerical stochastic perturbation theory.
In this presentation we show the results of the perturbative coefficients of the gradient flow coupling and its flow time dependence.
We investigate the beta function from the flow time dependence and
discuss the lattice artifacts in the large flow time in taking the large-$N$ limit.
