Assuming an underlying Hilbert space $\mathcal{H}$ of size $D$ for wavefunctions $|\psi\rangle$ of a many-body system, density matrices live the quadratically larger operator space $\hat{\rho}\in \mathcal{O}=\mathcal{H}\otimes\mathcal{H}^*$ of size $D^2$. This has strong limiting consequences, for example if one wants to solve master equations numerically. Of course, there are a few additional constraints (hermitian, positive definite, unit trace) to density matrices compared to general elements of $\hat{\rho}\in \mathcal{O}$; but I don't think they make too much conceptual difference.
Now, consider the purity of states $\mu=Tr \hat\rho^2 $ (or equivalently, any other notion of entropy). We know that pure states that can be described with a wavefunction alone, have $\mu=1$, a single scalar constraint. This means that the pure density matrices, live in a hypersurface of $\mathcal{O}$ with dimension $D^2-1$. This seems like a much less efficient description than the wavefunctions.
Now, consider the opposite. We start from the 'pure state' hypersurface and then leave it by tuning $\mu$. This could be done in a unique, orthogonal way (following the gradient of $\mu$) by using Gramm-Schmitt orthogonalization. This implies that you would only need one additional number to describe mixed states, as compared to pure.
Then, we could conclude that for all mixed states, we can describe them as wavefunctions in $D+1$ dimensions only! Actually, we could even drop the $+1$, by exploiting a remaining degree of freedom, the normalization of the wavefunction!
The trade-off of such construction would be that we lose the linearity of quantum mechanics, but still, the result, that mixed states can be described as wavefunctions, would seem spectacular to me.
Probably, I'm missing something in my reasoning, because otherwise this result would be quite known and exploited. Does anyone know what it is?