Solution (source code)

= Solution

Construct ensemble purifications
$$
|\Phi_1\rangle=\sum_j\sqrt{p_j}|\phi_j\rangle|j\rangle,
\qquad
|\Phi_2\rangle=\sum_k\sqrt{q_k}|\psi_k\rangle|k\rangle.
$$
They have the same reduced state exactly when $\rho_1=\rho_2$. By the <unitary freedom of purification>, this holds exactly when $|\Phi_2\rangle=(I\otimes U)|\Phi_1\rangle$ for a unitary $U$. Comparing reference-basis coefficients gives the <Hughston–Jozsa–Wootters theorem> relation
$$
\boxed{\sqrt{q_k}|\psi_k\rangle
=\sum_jU_{kj}\sqrt{p_j}|\phi_j\rangle}.
$$
Conversely, substituting this relation and using $U^\dagger U=I$ immediately gives $\rho_2=\rho_1$.