High-Order, High-Dimension, and Sparse Tensor (HOHDST) data originates from real industrial applications, i.e., social networks, recommender systems, bio-information, and traffic information. Sparse Tensor Decomposition (STD) can project the HOHDST data into low-rank space. In this work, a novel method for STD of Kruskal approximating the core tensor and stochastic strategy for approximating the whole gradient is proposed which comprises of the following two parts: (1) the matrization unfolding order of the Kruskal product for the core tensor follows the multiplication order of the factor matrix and then the proposed theorem can reduce the exponential computational overhead into linear one; (2) stochastic strategy adopts one-step random sampling set, the volume of which is much smaller than original one, to approximate the whole gradient. Meanwhile, this method can guarantee the convergence and save the memory overhead. Due to the compactness of the same order matrix multiplication and parallel access from stochastic strategy, the speed of cuFastTucker can be further reinforced by GPU. Furthermore, %because large-scale HOHDST data cannot be accommodated in a single GPU, a data division and communication strategy of cuFastTucker is proposed for data accommodation on Multi-GPU. cuFastTucker can achieve the fastest speed and keep the same accuracy and much lower memory overhead than the SOTA algorithms, e.g., P$-$Tucker, Vest, and SGD$\_$Tucker. The code and partial datasets are publically available on "https://github.com/ZixuanLi-China/FastTucker".