  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.
