As a high-order approximation model to Radiative Transfer Equation, simplify spherical harmonic (SPN) approximation has become a hot research topic in optical molecular imaging research. However, low computational efficiency imposes restrictions on its wide applications. This paper presented a graphics processing units (GPU)-parallel accelerated strategy for solving SPN model. The proposed strategy adopted compute unified device architecture (CUDA) parallel processing architecture introduced by NVIDIA Company to build parallel acceleration of two most time-consuming modules, generation of stiffness matrix and solving linear equations. Based on the feature of CUDA, the strategy optimized the parallel computing in tasks distribution, use of memory units and data preprocessing. Simulations on phantom and digital mouse model are designed to evaluate the accelerating effect by comparing the time for system matrix generation and average time of each step iteration. Experimental results show that the overall speedup ratio is around 30 times, which exhibit the advantage and potential of the proposed strategy in optical molecular imaging.