Suppose we are given information about distributions of random permutations $\sigma, \tau : \Omega \to S_n$ as follows:
$$p^1_{k,l} = \mathbb P(\sigma(k) = l), p^2_{k',l'} = \mathbb P(\tau(k) = l), p^{1,2}_{k,k'} = \mathbb P(\sigma(k) = \tau(k')),$$ for any $k,l,k',l' \leq n$.
We want to construct a "joint distribution", i.e. a tensor $Q \in \mathbb R^{n\times n \times n \times n}$ with the interpretation
$$Q_{k,l}^{k',l'} = \mathbb P(\sigma(k) = l, \tau(k') = l').$$
You may forget about probability now if you like. The actual marginal distributions $\mathbb P(\sigma = \pi)$ and $\mathbb P(\tau = \pi)$ are in fact not given and can be whatever you like as long as they satisfy the above equations.
In this case I'm just going to assume we are given doubly stochastic, symmetric matrices $p^1$, $p^2$, $p^{1,2}$ and we want to find $Q$, such that
\begin{gather*} \sum_{k=1}^n Q_{k,l}^{k',l'} = \sum_{l=1}^n Q_{k,l}^{k',l'} = p^2_{k',l'} \\ \sum_{k'=1}^n Q_{k,l}^{k',l'} = \sum_{l'=1}^n Q_{k,l}^{k',l'} = p^1_{k,l} \\ \sum_{l=1}^n Q_{k,l}^{k',l} = p^{1,2}_{k,k'} \end{gather*}
for all indices that we are not summing over.
At this point they are a lot more equations than unknowns. I feel like I should add a condition or two to make everything nicer. I just don't know what could be reasonable. In any case if I could decide I prefer "maximally dense" solutions over sparse solutions (I think this would make more sense for my application).
How can we find a "nice" solution this system of equations? I would think someone has studied and solved a problem like this before somewhere. I would like to know about a general methodology for approaching this kind of problem, because I want to extend it to multivariate distributions later if possible. Perhaps there is an abstract algebraic problem that generalizes this setting?
I'm not sure this question is quite up to standards, but I somehow always get completely overwhelmed trying to make sense of this problem, because the number of equations and unknowns I just so high and I really need to solve it for arbitrary $n$ and can't seem to find patterns for low $n$ solutions.