Hmm. I notice that the re-normalization trick offered by dpb sometimes seems to not converge well. That begs the question of whether there is a better algorithm. And, of course, there is. At least, there can be. The important question is if I can make it efficient for larger values of n.
A while ago when I last answered this question, I observed that you can represent all such doubly stochastic matrices as a convex linear combination of permutation matrices. (Of course, had I read the article given by dpb back then, it would have saved me some thought.)
The point is, that realization makes it easy to generate the set of such matrices, although for large n, that may become computationally intensive, at least by the initial algorithm I have in mind.
So, first, because we want here matrices with all zeros on the diagonal, that means we need to consider only permuation matrices with the same property. What do they correspond to? (Hint: this is called a derangement).
So let me see what happens for 3x3 matrices. Here are a couple of them.
0 0.528057118801456 0.471942881198544
0.471942881198544 0 0.528057118801456
0.528057118801456 0.471942881198544 0
0 0.700992660921403 0.299007339078597
0.299007339078597 0 0.700992660921403
0.700992660921403 0.299007339078597 0
You should seee a pattern in there. They are in fact a linear combination of the two matrices
a1 = [0 1 0;0 0 1;1 0 0]
0 1 0
0 0 1
1 0 0
a2 = [0 0 1;1 0 0;0 1 0]
0 0 1
1 0 0
0 1 0
Thus ALL 3x3 doubly stochastic matrices with zeros on the main diagonal can be written as a simple one dimensional family of matrices:
bistoc3 = @(t) t*a1 + (1-t)*a2;
So, just pick a random value of t in the interval [0,1],
0 0.3 0.7
0.7 0 0.3
0.3 0.7 0
The point being that for a 3x3 matrix, there need be no iterative scheme involved. No worries about convergence. And bistoc3 clearly will produce a family of matrices with uniformly distributed elements.
In a higher number of dimensions, this gets more complex, of course. The number of derangements of n numbers is given as the subfactorial of n, written as !n (as opposed to n! for the standard factorial.)
So let us look at the derangement idea, and what that does for us. For a 3x3 matrix, we need to look at derangements of 3 numbers. We can generate all of the derangements of the vector 1:3 as:
P3 = perms(1:3);
P3(any(P3 == 1:3,2),:) = 
3 1 2
2 3 1
Of course this works reasonably well, but 3 is a small number in context. How many derangements of 4 numbers are there? That is, how many 4x4 permutation matrices exist, with all zeros on the main diagonal?
P4 = perms(1:4);
P4(any(P4 == 1:4,2),:) = 
4 3 2 1
4 3 1 2
4 1 2 3
3 4 2 1
3 4 1 2
3 1 4 2
2 4 1 3
2 3 4 1
2 1 4 3
In fact, there are 9 derangements of the vector 1:4. This seems to be growing almost factorially. It should seem that way, because the number of derangements of a vector of length n is the subfactorial of n, written !n instead of the notation n! for the traditional factorial. It is not even that difficult to derive a recursive formula for the number of derangements, as
!n = (n-1)* (!(n-1) + !(n-2))
Of course we would start the recursion off with
!0 = 0, and !1 = 1.
These numbers grow rapidly, much like the factorial function itself, with !10 = 1334961.
The point of all this is we could write the set of random 4x4 doubly stochastic matrices with a zero main diagonal as a convex linear combination of the 9 corresponding derangement permutation matrices, as described in P4 above.
But this gets to the gist of the problem. For small n, say 3, 4, 5, it is not a problem. For example, how might I pose a fully uniform solution in the 4x4 case?
We can think of the permutation matrices as corner vertices of a polytope in a high number of dimensions.
There are 9 valid derangements of 1:4. This should work, as long as the dimension is not too large. (10x10 is probably a reasonable upper limit.)
eye4 = eye(4);
P4 = perms(1:4);
P4(any(P4 == 1:4,2),:) = ;
np4 = size(P4,1);
W = randfixedsum(np4,1,1,0,1);
A = zeros(4);
for i = 1:np4
A = A + W(i)*eye4(:,P4(i,:));
As you can see, this produces a valid matrix in the form you want. I used Roger Stafford's RANDFIXEDSUM code to compute a set of weights in 9 dimensions.
0 0.36499 0.35265 0.28236
0.29634 0 0.20877 0.49489
0.52519 0.25205 0 0.22276
0.17847 0.38295 0.43858 0
1 1 1 1
The question becomes what is the distribution of the elements? I have a funny feeling they are not as uniform as I want them to be, as the construction method I used fails because of the law of large numbers. Again, the elements of A look to be a bit to much close to the expected 1/3 average for my eyes.
So, can I write a better scheme? Yes, probably... (thinking on this.)