Solution (source code)

= Solution

The states $\rho_A^i=|a_i\rangle\langle a_i|$ are pure, so a purification is
$$
\boxed{
|\Psi\rangle
=\sum_{i=1}^n\sum_{\substack{j=1\\j\ne i}}^n
\sqrt{p_{ij}}\,
|i\rangle_{X_1}|j\rangle_{X_2}
|a_i\rangle_{A_1}|a_j\rangle_{A_2}
|ij\rangle_R}.
$$
Tracing out the orthonormal reference labels removes all cross terms and recovers the stated classical-quantum mixture. Its rank is the number of nonzero $p_{ij}$, so that is the minimum reference dimension for a particular state. The smallest dimension that can purify every state of the stated form is
$$
\boxed{\dim H_R=n(n-1)}.
$$