# How do I Generate Doubly-Stochastic Matrices Uniform Randomly?

A doubly-stochastic matrix is an $n\times n$ matrix $P$ such that

$\displaystyle\sum_{i=1}^n{p_{ij}}=1$

and

$\displaystyle\sum_{j=1}^n{p_{ij}}=1$

where $p_{ij}\ge 0$.

Can someone suggest an algorithm for generating these matrices uniform randomly?

-
Seen this? –  Guess who it is. Dec 1 '10 at 6:28
@J.M. Now, I have. (+1) –  Henry B. Dec 1 '10 at 7:26

What we want is to generate a bistochastic matrix according to the Haar measure, which is the unique distribution which is invariant to multiplication by bistochastic matrices from both sides.

The standard algorithm is to take some iid matrix (each entry is chosen iid from some distribution over non-negative numbers) and then repeatedly make it row-stochastic and column-stochastic - this is like projecting the matrix to the linear subspace of stochastic matrices. The process converges pretty fast, and we get some random bistochastic matrix, even though not according to the wanted distribution.

A recent paper discusses these issues in the context of estimating the measure of all bistochastic matrices (an open problem).

-
Nice. I just finished implementing it. Thanks! –  Henry B. Dec 1 '10 at 7:26

You can first generate a random unitary matrix and then square the absolute values of all entries. Here is Mathematica code that does this:

(* Random real and complex numbers with normal distribution *)
RR := RandomReal[NormalDistribution[0, 1]];
RC := RR + I*RR;
(* Random matrix from Ginibre ensemble *)
RG[n_] := Table[RC, {n}, {n}];
(* Random unitary matrix *)
RU[n_] := Orthogonalize[RG[n]];
(* Random doubly stochastic matrix *)
RDS[n_] := Abs[RU[n]]^2;


Run RDS[5] to generate a $5 \times 5$ random doubly stochastic matrix.

-

Suppose you want a random rational point in the Birkhoff polytope (set of double stochastic matrices), with denominators dividing $N$. This is equivalent to a random contingency table with all row and column sums equal to $N$. A reversible Markov chain exists on the set of tables and is rapidly mixing. This can be used for approximate sampling by following the chain for a large number of steps. The limit distribution of the chain is close to uniform because away from the boundary of the polytope, all vertices of the chain have the same degree. Exact sampling from the limit distribution of the chain is possible, but because the limit distribution is not quite uniform this is not the same as obtaining an exact sample from the uniform distribution on contingency tables. The latter is a hard unsolved problem as far as I know.

-
Can you add a link to the rapidly mixing Markov chain? (and possibly also explain how it goes) –  Yuval Filmus Dec 1 '10 at 6:41
The chain is generated by moves where to a 2x2 submatrix you add or subtract the 2x2 matrix [(+1, -1) (-1,+1)] or its negative. I don't know the reference for rapid mixing but am pretty sure you can find it (for general contingency tables) in papers of Diaconis and his collaborators. –  T.. Dec 1 '10 at 7:21