The case in which the initial state is bipartite and pure has a geometric solution, please see the original construction by Życzkowski and Kuś and the review: by Bengtsson and Życzkowski. It also appears in section 16.4 of the book: GEOMETRY OF QUANTUM STATES An Introduction to Quantum Entanglement.
The characterization of the bipartite pure states having the same marginals is that they lie on the same fiber in the bundle structure of a space of equally entangled states over the space of unitary equivalent density matrices whose nonvanishing eigenvalues are the squares of the original state Schmidt coefficients. This characterization can be explained as follows:
For a bipartite system, the non-vanishing eigenvalues of the reduced density matrix are the squares of the non-vanishing Schmidt coefficients of the initial pure states. The action of a local unitary transformation $U \otimes V$ on the initial state produces a unitarily displaced reduced density matrix by $U$ on $\mathcal{H}_1$ and $V$ on $\mathcal{H}_2$. Thus, the space unitarily equivalent reduced density matrices with a given set of non-vanishing eigenvalues is obtained by tracing the states belonging to the space of equally entangled bipartite pure states, whose nonvanishing Schmidt coefficients are the square roots of the eigenvalues. By equally entangled states, it is meant that they lie on a local unitary transformation orbit.
Said differently, a space of equally entangled bipartite pure states is a fiber bundle over the unitary orbit of the reduced density matrix whose non-vanishing eigenvalues is equal to the initial state Schmidt vector.
The fibers are isomorphic to $U(N)/U(m_0)\times U(1)$, where $N$ is the dimension of the Hilbert space which was traced on and $m_0$ is the number of vanishing Schmidt coefficients. The fibers corresponds to spaces of states having the same reduced density matrix. Please see table 16.1 in the third reference for the various orbit structures in the case of 2 qubits and 2 qutrits.
Update
Given a pure state $|\Psi\rangle$ living in $\mathcal{H}_1\otimes\mathcal{H}_2=\mathbb{C}^n\otimes\mathbb{C}^n $, then its local unitary orbit is
$$\mathcal{O} = \mathrm{Span} \left(U(n)\otimes U(n) \circ \Psi \right)$$
This is the total space of the bundle: $\mathcal{F}\rightarrow \mathcal{O}\rightarrow \mathcal{M}$.
Suppose that the Schmidt coefficients of this state are given by:
$\mathbf{v_s} = \{\underbrace{0,…}_{m_0 \mathrm{times}}, \underbrace{\lambda_1,…}_{m_1\mathrm{times}}, \underbrace{\lambda_2,…}_{m_2\mathrm{times}}, …\}$
Then (after tracing on the second factor) base space of the bundle is:
$$\mathcal{M} = \mathrm{Span} \left(U(n) \circ \mathrm{tr}_2\left(|\Psi\rangle\langle|\Psi|\right) \right)= \mathrm{Ad}\left(U(n)\right) \mathrm{diag}(|\mathbf{v_s}|^2)= U(n)/U(m_0)\times U(m_1)\times U(m_2)…$$
The fiber of the bundle is:
$$\mathcal{F} = U(n)/U(m_0)$$
(The reason that $U(m_0)$ appears in the isotropy group of both the base space and the fiber is that the isotropy group of a unit matrix $\mathbf{1}$ is $U(m_1)$: $\mathbf{1} = U \mathbf{1}U^{\dagger} $, while the isotropy group of a zero matrix is $U(m_0) \times U(m_0)$: $\mathbf{0} = U \mathbf{0}V^{\dagger} $)
The span in the base space is due to the first $U(n)$ factor of the local unitary group $ U(n)\otimes U(n)$, i.e., the factor acting on $\mathcal{H}_1$, while the span in the fiber is the second $U(n)$ factor acting on $\mathcal{H}_2$. The base space $M$ is the unitary orbit of the reduced state, thus, if we wish to characterize states descending to the same reduced state, we have to fix the second $U(n)$ factor, thus remain with a space isomorphic to the fiber $\mathcal{F}$.