A slightly different spin on the characterisation of the set of ensemble decompositions of a given state.
Let $\rho$ be a generic positive semidefinite operator. I'm not going to assume normalisation here, as it's not necessary for the argument. Being it positive semidefinite, $\rho$ admits an eigendecomposition
$$\rho=\sum_k p_k \mathbb{P}(\psi_k),$$
where $p_k>0$ and $\mathbb{P}(\psi_k)\equiv|\psi_k\rangle\!\langle\psi_k|$ are the rank-1 projections on the states $|\psi_k\rangle$.
Note in particular that by definition of eigendecomposition we have $\langle\psi_j,\psi_k\rangle=\delta_{jk}$.
In general we can have $\sum_k p_k\neq1$ if $\rho$ is not normalised. If $\operatorname{tr}(\rho)=1$ then $\sum_k p_k=1$ and this equation expresses $\rho$ as a convex combination of rank-1 operators.
More generally, this equation expresses $\rho$ as an affine combination of such operators. In the following, just replace all mentions of "affine" with "convex" to specialise to the $\operatorname{tr}(\rho)=1$ case.
The goal is to characterise all possible ways to decompose $\rho$ as $\rho=\sum_k \lambda_k \mathbb{P}(\phi_k)$ for some set of $\lambda_k>0$ and not-necessarily-orthogonal states $|\phi_k\rangle$.
The trick is to observe that the non-negativity of the coefficients allows to write $\rho=\sum_k \mathbb{P}(\sqrt{\lambda_k}\phi_k)$, which can in turn be written concisely as $\rho= AA^\dagger$ with
$$A \equiv \sum_k \sqrt{\lambda_k} |\phi_k\rangle\!\langle k|.$$
In words, $A$ is a linear operator whose columns are the vectors $\sqrt{\lambda_k}|\phi_k\rangle$.
Let me stress that not only any decomposition of $\rho$ in terms of rank-1 operators has this form, but also any linear operator $A$ such that $\rho=AA^\dagger$ corresponds to one such decomposition.
In other words, the problem we're discussing is mathematically equivalent to that of characterising the set of linear operators $A$ such that $\rho=AA^\dagger$, for a given $\rho\ge0$.
Note that we also don't need $A$ to have the same dimension of $\rho$.
If $\rho$ is $d$-dimensional, meaning $\rho\in\operatorname{Lin}(\mathbb{C}^d)$, we can have $A$ of the form $A\in\operatorname{Lin}(\mathbb{C}^n,\mathbb{C}^d)$ for any $n\ge \operatorname{rank}(\rho)$. The idea is here that the dimension of the domain of $A$ corresponds to the number of terms in the decomposition, which can't be smaller than the rank of $\rho$, but can be arbitrarily large (in fact, it can infinite, though I'm only considering finite-dimensional cases here).
So how do we find the linear operators $A$ such that $\rho=AA^\dagger$? The answer is easy thinking in terms of SVDs.
This relation imposes $A$ to have singular values equal $\sqrt{p_k}$, and left principal components equal to the eigenvectors of $\rho$.
Or more explicitly, the operators we are looking for are all and only those of the form
$$A = \sum_k \sqrt{p_k} |\psi_k\rangle\!\langle u_k|,$$
where $|\psi_k\rangle$ are again the eigenvectors of $\rho$. This sum thus clearly contains a number of elements equal to the rank of $\rho$.
Furthermore, $\{|u_k\rangle\}$ is here an orthonormal basis for the domain space of $A$, and any such orthonormal basis gives a valid $A$. And any decomposition of $\rho$ corresponds to such a choice of basis.
This freedom can be expressed concisely writing that
$$A= \sqrt\rho V^\dagger,$$
where we defined the isometry $V\equiv \sum_k |u_k\rangle\!\langle \psi_k|$. Note that this is an isometry, and not in general a unitary operator, because $|u_k\rangle$ can live in a space of higher dimension than $|\psi_k\rangle$.
In conclusion, the set of decompositions for $\rho$ corresponds to the set of linear operators of the form $A=\sqrt\rho V^\dagger$ where $V$ spans all possible isometries. The elements in the decomposition of $\rho$ correspond to the columns of $A$.
More explicitly this means that $\rho$ decomposes as $\rho=\sum_k a_k a_k^\dagger$ where
$$a_k \equiv A|k\rangle = \sum_j \sqrt{p_j}|\psi_j\rangle \overline{\langle k|V|\psi_j\rangle}
= \sum_j \sqrt{p_j} |\psi_j\rangle \bar u_{kj},$$
defining for convenience $U$ as the isometry such that $\langle k|U|j\rangle=\langle k|V|\psi_j\rangle$.
This is the formula one usually finds written explicitly.