Random permutations by probability matrix

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.

• Welcome to Computer Science! Note that you can use LaTeX here to typeset mathematics in a more readable way. See here for a short introduction. – FrankW Oct 12 '14 at 19:29
• What are the permutation probabilities for the matrix $$\begin{matrix}1/3&1/3&1/3\\1/2&1/2&0\\1/2&1/2&0\end{matrix}$$? – David Eisenstat Oct 12 '14 at 20:17
• For this matrix are possible only following permutations: 312, 321, both with same frequency. – michal Oct 13 '14 at 10:52
• A mixture of 312 and 321 would correspond to the matrix $\begin{bmatrix} 0&0&1\\1/2&1/2&0\\1/2&1/2&0 \end{bmatrix}$. – Yuval Filmus Oct 13 '14 at 13:22
• @YuvalFilmus Another way to interpret the matrix entries would be to sample each permutation $\pi$ with probability proportional to $\prod_{i=1}^n p_{i,\pi(i)}$. I'm not sure whether this interpretation differs in the bistochastic case. – David Eisenstat Oct 13 '14 at 15:43

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.
• 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 produce all six possible permutations with same frequency. But my problem is opposite. I want to generate random permutations by predefined probability matrix, not create probability matrix from permutations. – michal Oct 13 '14 at 9:42
• @michal No, you misunderstood me completely. I can't understand your model. In my model, $p_{ij}$ is the probability that $i$ goes to $j$. Under this model, you can generate only and all bistochastic matrices. What happens in your model, I can't say, since I don't understand it. Both examples you give are bistochastic, so I suspect our models are actually the same. – Yuval Filmus Oct 13 '14 at 13:00