Suppose that $X$ and $Y$ are finite-dimensional Hilbert spaces with $\mathrm{dim}(X)=n$ and $\mathrm{dim}(Y)=m$. Let $\rho$ be a density operator on $X \otimes Y$ such that the reduced density operator on $X$ is a classical state; that is, $\mathrm{Tr}_{Y}(\rho) = \sum_{i=1}^n p(i) u_i u_i^*$, where $\{u_i\}$ is an orthonormal basis of $X$, $\sum_{i=1}^n p(i) = 1$, $p(i) \geq 0$ for all $i$, and $u_i^*$ is complex conjugate of $u_i$.
How can we prove that the compound state $\rho$ is always of the following form: $\rho = \sum_{i=1}^n p(i) u_i u_i^* \otimes \sigma_i$, where $\sigma_i$ is an arbitrary density operator on $Y$ for all $i$?