Random permutations by probability matrix

Question

I have the following problem:

I need to generate $\ell$ random permutations each of length $n$ from a list of $m$ elements ($m \ge n$) by a predefined probability matrix $P$ of size $n$ x $m$.

Distribution probability matrix $P$ = {$p_{i,j}$} describes the probability that i-th position is occupied by $j$-th element ($i\in \{1,\dots,n\}$ and $j \in \{1,\dots,m\}$).

Of course $1 \ge p_{i,j} \ge 0$ and $\sum\limits_{j=1}^m(p_{i,j}) = 1$ (rows of matrix P are normalized, to be equal 1)

$\sum\limits_{i=1}^n(p_{i,j}) > 0$ (each element has nonzero probability to be at least on one position)

If $p_{i,j} = 0$, then $i$th position of permutation can not be occupied by $j$th element. If $p_{i,j} = 1$, then $i$th position is occupied only by $j$th element

Example: For $$P = \begin{bmatrix} 0.5 & 0.5 & 0 \\ 0 & 0.5 & 0.5 \\ 0.5 & 0 & 0.5\\ \end{bmatrix},$$ I get only two permutations [1 2 3] and [2 3 1] with the same frequency. But for $$P = \begin{bmatrix} 0.01 & 0.99 & 0 \\ 0 & 0.01 & 0.99 \\ 0.99 & 0 & 0.01\\ \end{bmatrix},$$ I significantly prefer permutation [2 3 1]

Realistic values $\ell = 1\,000$–$10\,000$, $n = 30$, $m = 50$

Is there any suitable and effective (fast) algorithm?

Edit: The problem is motivated by permutation sampling for stochastic combinatorial optimization. Permutations are parametrized via probability matrix P. The probability matrix P is generated by black-box objective function responses.

Edit2: This is typical algorithm for permutation sampling at stochastic optimization problem:

Generation of partial permutations $x$ by probability matrix $P$

1. Generate a random permutation $(\pi_1, ... ;\pi_{n})$ of the set of positions $(1,...,m)$.

2. Define $P(1) = P$ and set $a = 1$.

3. Generate $x_{\pi_{a}}$ according to the distribution formed by the $a$-th row of $P(a)$ , that is $(p_{(\pi_{a},1)},... ,p_{(\pi_{a},m)})$. Thus, element $x_a$ is placed into position $a$.

4. Obtain $P(a+1)$ from $P(a)$ by first setting the $x_{a}$-th column of $P(a)$ to 0, and then renormalizing the rows to sum up to 1.

5. If $a = n$ then stop and return $x$; otherwise set $a = a + 1$ and return to step 3.

Matrix $P$ is during optimization cycle obtained by update algorithm, which produce stochastic (not bistochastic) matrix. I want to be sure, that this algorithm is mathematically correct.

Answer 1

Your problem is not well-defined. As David Eisenstat notes in his comment, your matrix actually has to be bistochastic rather than just stochastic, since every convex combination of permutation matrices is bistochastic. Also, there could be many ways of representing a given bistochastic matrix as a convex combination of permutation matrices. As a simple example, the matrix $$\begin{bmatrix} 1/3 & 1/3 & 1/3 \\ 1/3 & 1/3 & 1/3 \\ 1/3 & 1/3 & 1/3 \end{bmatrix}$$ can be obtained by taking all permutations or just $123,231,312$.

The Birkhoff–von Neumann theorem states that every bistochastic matrix is a convex combination of permutation matrices. Its algorithmic proof gives a polynomial time algorithm that gives such a combination with at most $n^2$ matrices, which is probably optimal or close to optimal. In the literature you might be able to find faster algorithms; if you do, please let us know.