The population of space debris grows to over a hundred million objects, making their orbital determination computationally expensive. This work proposes an efficient propagation method for space debris. The initial orbital states of massive debris within a bounded region are characterized by some special samples, which are propagated to any future time. Then the orbital states of debris are constructed by Chebyshev polynomial interpolation at the propagated samples. A sparse sampling strategy is proposed to enable the sample number to be much smaller than the space debris, and a parallel integrator is implemented to guarantee efficient propagation. The performance of the proposed method is compared to previous propagation methods for the propagation of debris in LEO with J2 perturbation. The effect on the accuracy and efficiency of polynomial order is discussed.