How large is the space of density matrices, really? 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?
 A: I'll try to formulate a physicist answer with a limited mathematical knowledge. First of all, density matrices (unlike wavefunctions) do not live in a vector space, but in the convex set $\mathcal{S}$ defined by:
$$
\mathcal{S}=\{\rho\in\mathcal{T}(\mathcal{H)}\text{ such that }\rho\geq 0,\text{Tr}[\rho]=1\},
$$
where $\mathcal{T}(\mathcal{H})$ is the space of operators on the Hilbert space $\mathcal{H}$ with well-defined trace. Therefore, I don't know if the "dimension" of the space of the density matrix can be rigorously defined. I see that there are some definitions of dimension of convex sets, but I can say that I have never seen them applied in the context of density matrices. This does not mean that mathematical physicists do not employ this notion, but anyway I think the issue at point here is the number of real parameters we need to define a density matrix. We can use the concept of generalized Bloch sphere [1]: any density matrix $\rho$ of a $N$-dimensional Hilbert space can be written as:
$$
\rho=\frac{\mathbb{I}}{N}+\sum_{i=1}^{N^2-1}\tau_i\sigma_i,
$$
where $\sigma_i$ are the generators of $SU(N)$ and $\{\tau_i\}_{i=1}^{N^2-1}$ is a collection of $N^2-1$ real numbers. That is to say, any density matrix is uniquely defined by $N^2-1$ real parameters. Some constraints (inequalities) must then be applied to this set of real numbers in order for $\rho$ to be a well-defined density matrix (see for instance Ref.$~$[2]), but this does not reduce the number of free real parameters. This seems to indicate that it is impossible to describe a density matrix through $2N$ real parameters only as you are suggesting.
What may be wrong in your reasoning? I'd say there is some trick about mathematical dimensions, which is not the dimension of a vector space anymore. Anyway, I didn't understand how to obtain the construction of any mixed state you are proposing, based on purity and a single pure state only.

This construction is clear only for the case $N=2$: in this scenario the set of physical states is given by the closed ball defined by $||\boldsymbol{\tau}||=\sqrt{\sum_{i=1}^3 \tau_i^2}<1$ (this is the only additional constraint on the collection $\{\tau_i\}_{i=1}^3$), and the set of pure states coincides with the surface of such ball. The purity of $\rho$ is a function of the length of the vector $||\boldsymbol{\tau}||$, $\mu[\rho]=\frac{1}{2}\left(1+||\boldsymbol{\tau}||^2\right)$. Therefore, any mixed state living inside the ball can be uniquely defined by its purity and the closest pure state on the surface. For $N>2$, however, many issues arise: even if the purity of any state $\rho$ is still given by the formula $\mu[\rho]=\frac{1}{2}\left(1+||\boldsymbol{\tau}||^2\right)$, the space of the physical states is a not a hyperball anymore and the set of pure states does not coincide with the boundary of the generalized Bloch sphere [2].
[1] Bengtsson and Zyczkowski, Geometry of quantum states: an introduction to quantum entanglement. Cambridge University Press (2006).
[2] Goyal et al., J. Phys. A: Math. Theor. 49, 165203 (2016).
A: If I understand the question, it seems to be about trying to find a hilbert space whose (pure) states encode the same information as (possibly mixed) density matrices of some other physical system.
I'd say this doesn't seem promising, for the following geometric/topological reason. Take the space of rays in the hilbert space (i.e. states up to scale). For an $n$-dimensional Hilbert space, that is the manifold
$$
\mathbb{CP}^{n-1}
$$
(check: for $n=2$ you get $\mathbb{CP}^1=S^2$, the Bloch sphere). That is a compact manifold (basically: finite volume), without boundary.
On the other hand, the space of density matrices is always a manifold with boundary! (This is basically some kind of filled ball. The boundary is given by the equation ${\rm Tr}[\rho^2]=1$.) For $n=2$, this boundary is just the Bloch sphere.
So while you can obviously pick the dimensionalities to make the manifold of mixed states have the same dimension as the manifold of pure states of some other system (that is essentially the counting you're doing in your question), the manifolds will always be different. And you should be concerned by this: if you compare the manifold of rays in hilbert space to the manifold of pure density matrices instead, you get the same manifolds on both sides. (In fact that shows how for pure states the density matrix and ket descriptions are equivalent.)
A: If $d$ is the dimension of your Hilbert space then the space of density matrices will have dimension $d^2-1$.
If $\{\vert\alpha\rangle,\alpha=1,\ldots, d\}$ span your Hilbert space, then you can make $d^2$ operators of the type $\vert \alpha \rangle\langle\beta\vert$, and you can expand the density matrix in terms of these operators.  The expansion coefficients can be complex i.e. $2d^2$ real numbers, but the density matrix must be hermitian:  the $d$ diagonal entries must be real, so that’s $d$ constraints,  and hermiticity bring another $\frac{1}{2}d(d-1)$ constraints on the off-diagonal complex numbers,or $d(d-1)$ constraints on real numbers.  Thus we have
$$
2d^2-2\times \frac{1}{2}d(d-1)-d= 2d^2-d^2+d-d=d^2\, .
$$
The last constraint is the normalization since the entries on the diagonal must sum to $1$.
The number of parameters grows quadratically, which is still polynomial, but somewhat expensive if you have large Hilbert spaces arising from multi-particle states.  This growth is one reason why Wigner functions are of interest as the number of parameters in the WF does not depend on the number of particles.
A: The density matrix is a function of the wave function. Assuming D is the complex dimension, given the freedom of multiplying by a phase $e^{i\alpha}$ for every dimension in the ensemble of states for the same density matrix, it should just be $D-1$.
I seem to have misunderstood OP's question. If it was for the calculation of possible density operators for unknown systems with a given Hilbert space, the number of such operators is the same as the number of "frame functions", which is defined on the Bloch sphere subject to $\sum_i f(x_i)=1$ for any orthonormal basis $x_i$, so it is infinite dimensional. You can see the description of frame function in Gleason's theorem :https://en.wikipedia.org/wiki/Gleason%27s_theorem.
