Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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




where $p_{ij}\ge 0$.

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

share|cite|improve this question
Seen this? – J. M. Dec 1 '10 at 6:28
@J.M. Now, I have. (+1) – Henry B. Dec 1 '10 at 7:26
up vote 7 down vote accepted

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

share|cite|improve this answer
Nice. I just finished implementing it. Thanks! – Henry B. Dec 1 '10 at 7:26

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.

share|cite|improve this answer
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

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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.