One of the most remarkable basis of the gravity data inversion is the recognition of sharp boundaries between an ore body and its host rocks during the interpretation step. Therefore, in this work, it is attempted to develop an inversion approach to determine a 3D density distribution that produces a given gravity anomaly. The subsurface model consists of a 3D rectangular prisms of known sizes and positions and unknown density contrasts that are required to be estimated. The proposed inversion scheme incorporates the Cauchy norm as a model norm that imposes sparseness and the depth weighting of the solution. A physical-bound constraint is enforced using a generic transformation of the model parameters. The inverse problem is posed in the data space, leading to a smaller dimensional linear system of equations to be solvedand a reduction in the computation time. For more efficiency, the low-dimensional linear system of equations is solved using a fast iterative method such as Lanczos Bidiagonalization. The tests carried out on the synthetic data show that the sparse data-space inversion produces blocky and focused solutions. The results obtained for the 3D inversion of the field gravity data (Mobrun gravity data) indicate that the sparse data-space inversion could produce the density models consistent with the true structures.