How do I Generate Doubly-Stochastic Matrices Uniform Randomly?

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).

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.

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.