A system of linear equations with way too many unknowns — constructing a bivariate distribution from marginals and "the diagonal" 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.
 A: $\newcommand\si\sigma$You are considering wrong unknowns, that is, a wrong "joint distribution".
The correct joint distribution of $\si$ and $\tau$ is given by their joint probability mass function (pmf) $p$ defined by the formula
$$p(s,t):=P(\si=s,\tau=t)$$
for $(s,t)\in S_n\times S_n$. The pmf $p$ must satisfy the restrictions
\begin{align}
\sum_{(s,t)\in S_n\times S_n} p(s,t)\,1(s(k)=l)=p^1_{k,l}, \tag{1}\label{1} \\ 
\sum_{(s,t)\in S_n\times S_n} p(s,t)\,1(t(k)=l)=p^2_{k,l}, \tag{2}\label{2} \\ 
\sum_{(s,t)\in S_n\times S_n} p(s,t)\,1(s(k)=t(l))=p^{1,2}_{k,l} 
\tag{3}\label{3}
\end{align}
for all $k,l$ in $[n]:=\{1,\dots,n\}$. Let $P:=P_{p^1,p^2,p^{1,2}}$ denote the set of all the joint pmf's $p$ satisfying these restrictions. One may note here is that the system \eqref{1}--\eqref{3} is a system of $3n^2$ linear (but not linearly independent) equations with $(n!)^2$ unknowns.
If $P$ is nonempty, you may want to select a joint pmf $p_{\min}\in P$ which (quasi-)minimizes whatever sparsity characteristic of $p\in P$ you had in mind.
Finally, if you wish, you can compute
$$Q_{k,l}^{k',l'}=P(\si(k)=l,\tau(k')=l')
=\sum_{(s,t)\in S_n\times S_n} p_{\min}(s,t)\,1(s(k)=l,t(k')=l')$$
for $k,l,k',l'$ in $[n]$.

One should not write things like $\tau\in S_n$. Indeed, a random permutation $\tau$ is, not a permutation, but a map into the set $S_n$ of permutations.
