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 [1].
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.