In the case of perfect spherical symmetry we can learn nothing about the density profile (other than that it is spherically symmetric). This is a consequence of Newton's shell theorem (or Birkhoff's theorem in GR) that says the potential (or metric) outside all spherically symmetric, static mass distributions of a given total mass are the same and spherically symmetric.
For all other cases, which includes every real object in the universe, there is no exact spherical symmetry and so, for example, an orbiting probe can yield detailed information about the distribution of mass inside the object it is orbiting. This is where much of the information we have about planetary interiors in our Solar System comes from.
Data from orbiting spacecraft can be used to constrain the mass multipole moments - the gravitational potential can be expressed in terms of a series expansion: mass monopole, (no dipole, because no negative masses), quadrupole etc. These in combination with physical modelling constrain the interior structure and density profile (e.g. Mazarico et a. 2014; Ni 2020).
Gravitational waves also contain information about the radial density profile of the (non-black hole) objects in a merging binary system. When the components get close they are distorted by tidal fields and this distortion, which depends on the interior equation of state and density profile, alters the gravitational wave emission and so can be used to model, for example, neutron star interiors (Chatziioannou et al. 2015.
I doubt there are uniquely invertible solutions, even if one were to measure all the multipoles to arbitrarily higher orders with high precision. It is usually the case that other physcal constraints are applied (e.g., hydrostatic equlibrium and an equation of state).