We have evaluated perturbation coefficients of Wilson loops up to $O(g^8)$ for the four-dimensional twisted Eguchi-Kawai model using the numerical stochastic perturbation theory (NSPT) in arXiv:1902.09847. In this talk we present a progress report on the higher order calculation up to $O(g^{63})$, for which we apply a fast Fourier transformation (FFT) based convolution algorithm to the multiplication of polynomial matrices in the NSPT aiming for higher order calculation. We compare two implementations with the CPU-only version and the GPU version of the FFT based convolution algorithm, and find a factor 9 improvement on the computational speed of the NSPT algorithm with SU($N=225$) at $O(g^{31})$. The perturbation order dependence of the computational time, we investigate it up to $O(g^{63})$, shows a mild scaling behavior on the truncation order.