Gravity inversion allows us to constrain the interior mass distribution of a planetary body using the observed shape, rotation, and gravity. Traditionally, techniques developed for gravity inversion can be divided into Monte Carlo methods, matrix inversion methods, and spectral methods. Here we employ both matrix inversion and Monte Carlo in order to explore the space of exact solutions, in a method which is particularly suited for arbitrary shape bodies. We expand the mass density function using orthogonal polynomials, and map the contribution of each term to the global gravitational field generated. This map is linear in the density terms, and can be pseudo-inverted in the under-determined regime using QR decomposition, to obtain a basis of the affine space of exact interior structure solutions. As the interior structure solutions are degenerate, assumptions have to be made in order to control their properties, and these assumptions can be transformed into scalar functions and used to explore the solutions space using Monte Carlo techniques. Sample applications show that the range of solutions tend to converge towards the nominal one as long as the generic assumptions made are correct, even in the presence of moderate noise. We present the underlying mathematical formalism and an analysis of how to impose specific features on the global solution, including uniform solutions, gradients, and layered models. Analytical formulas for the computation of the relevant quantities when the shape is represented using several common methods are included in the Appendix.