Low calculation efficiency and insufficient resolution in depth direction exist in inversion of underground 3D density distribution. In this paper, we proposed a fast inversion algorithm. It decomposed the gravity anomaly on multi-scale with wavelet, represented the original data sparsely with wavelet coefficients on each scale and carried on the inversion in wavelet domain. An appropriate threshold is set to process the coefficients in order to enhance the sparsity of coefficient matrix. This would furthermore improve the compressed ratio of data and save calculation time. Gravity anomalies on each scale represent response of sources at different depth as inverse relation exists between scale and frequency. So the inversion results would mainly reflect density distribution at the corresponding depth and the final inversion could be achieved by summing up results at all scales. It is not necessary to set the range of depth related to anomalies at each scale and besides, the inversion scheme is applied without depth weighting. The method could increase the resolution in depth direction of inversion results effectively and provide more detailed deep density distribution. As an iterative algorithm which can take advantage of sparse matrix to improve calculation efficiency, conjugate gradient was used for inversion. The proposed method will be applied into inversion of synthetic model data and real gravity anomaly to demonstrate its effect.