Can the wavefunction of a system be reconstructed from density matrices of subsystems? Suppose we have several interacting particles in pure state $\left|\psi\right>$. For each of particles we can extract density matrices via
$$\rho_i(x_i,x_i^\prime)=\int \left<X_i,x_i\middle|\psi\right>\left<\psi\middle|X_i,x_i^\prime\right>dX_i,$$
where $x_i$ is coordinates of particle $i$ and $X_i$ is coordinates of the rest of system.
Can the original wave function $\left<X\middle|\psi\right>$ of the whole system be restored from all the $\rho_i$?
 A: No, it is not possible in general, because the particles are almost surely to be entangled (since they are interacting, as mentioned).
The reason is that each $\rho_i$ is nothing but the reduced density operator of the $i$ th particle. Indeed, in the form you have written, 
$$\rho_i(x_i,x_i') = \langle x_i | \rho_i | x_i' \rangle = \langle x_i | \mathsf{Tr}_{X_i \neq x_i}(|\psi \rangle \langle \psi |) | x_i' \rangle ,$$
where the $X_i$ go over the rest of the system.
Even if all the $\rho_i$ are known, one cannot find $\rho \equiv |\psi \rangle \langle \psi |$, the state of the system, because there can exist correlations among the variables of different particles that cannot permit one to write the following:
$$
\rho = \rho_1 \otimes \rho_2 \otimes .. \rho_i \otimes .. \ \ \ \ \mathsf{(Incorrect)}
$$
This is just the continuous variable version of entanglement. As an example with qubits, consider the Bell state 
$$
|\Psi \rangle = \frac{|01\rangle + |10\rangle}{\sqrt{2}}
$$
and trace out one of the qubits. This yields
$$
\rho_1 = \rho_2 = \frac{|0\rangle\langle 0|+|1\rangle\langle 1|}{2},
$$
But even a state with only classical correlations, 
$$
\rho = \frac{|00\rangle \langle00| +|11\rangle\langle11|}{2}
$$
would give the same reduced density operators, this state being much different from the entangled one considered. This shows that there is more information in the full state $\rho$ than in the partial knowledge of each of the subsystems $\rho_i$.
A: As Abhinav states, in general, it is not possible. Assuming low interaction you can however start with a Weyl transform
$$\begin{align*}
  P_i(x,p) &= \frac{1}{\pi\hbar}\int_{-\infty}^\infty \langle x+y| \hat{\rho}_i |x-y \rangle e^{-2ipy/\hbar}\,dy &\Bigg|\quad \hat\rho_i =\iint du\,dv\,\rho_i(u,v)|u\rangle\langle v|
\\ &= \frac1{\pi\hbar}\int_{-\infty}^\infty \rho_i(x+y,x-y) e^{-2ipy/\hbar}\,dy
\end{align*}$$
Then use $|\psi_i(x)|^2 = \int_{-\infty}^\infty P_i(x,p)\,dp$ and $|\phi_i(p)|^2=\int_{-\infty}^\infty P_i(x,y)\,dx$ to obtain the intensity of the individual wave functions in coordinate- and momentum-space. Finally, reconstruct the phase e.g. via Gerchberg-Saxton. Obviously this will fail horribly once sufficient interaction/mixture occured...
