Generate random Bistochastic matrix with zero diagonals

Samuel Spedding (view profile)

on 2 Sep 2019
Latest activity Edited by dpb

on 3 Sep 2019

dpb (view profile)

I'm looking to generate a load of n x n random matrices with non-negative components and with zeros on the diagonals, and whose rows and columns all sum to 1. So far I've managed to create matrices whose rows sum to 1, but I'm struggling with how to get the columns to sum to 1 also.
F=zeros(n);
for i = 1:n
for j = 1:n
if i~=j
F(i,j) = exp(randn);
end
end
F(i,:) = F(i,:)/sum(F(i,:));
end
Any help would be gladly appreciated. Thanks.

R2017b

on 3 Sep 2019
Edited by dpb

dpb (view profile)

on 3 Sep 2019

function a=bistochastic(a,tol)
% return bistochastic matrix from input nonnegative matrix
if nargin<2, tol=0.001; end
s1=sum(a,2);
s2=sum(a);
while ((abs(max(s1)-1))>tol) | ((abs(max(s2)-1))>tol)
a=a./s1;
s1=sum(a,2);
s2=sum(a);
a=a./s2;
end
end
will converge if begin with matrix with arbitrary nonnegative entries.
Example:
>> A=rand(4).*not(eye(4)); % initial random array w/ zero diagonal for input
>> A
A =
0 0.4505 0.1524 0.0782
0.6541 0 0.8258 0.4427
0.6892 0.2290 0 0.1067
0.7482 0.9133 0.9961 0
>> [sum(A); sum(A,2).']
ans =
2.0914 1.5929 1.9743 0.6275
0.6811 1.9226 1.0248 2.6576
>> B=bistochastic(A,1E-5); % create bistochastic from it w/ 10^-5 error in sum
>> B =
0 0.5038 0.2195 0.2768
0.2062 0 0.3425 0.4513
0.5436 0.1844 0 0.2720
0.2501 0.3118 0.4380 0
>> [sum(B); sum(B,2).']
ans =
1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000
>>
You could, of course, have the function generate the initial matrix; I left it general so pass in the special condition of the zero diagonal. Alternately, that could be an input argument instead.
I ran across the above years ago...see <Random bistochasitic matrices> for some more.

John D'Errico

John D'Errico (view profile)

on 3 Sep 2019
It is actually an interesting problem. I hacked your code slightly.
function a=bistochastic(n,tol)
a = rand(n);
a = a - diag(diag(a));
s1=sum(a,2);
s2=sum(a);
iter = 0;
itmax = 100000;
while ((abs(max(s1)-1))>tol) || ((abs(max(s2)-1))>tol)
iter = iter + 1;
a=a./s1;
s1=sum(a,2);
s2=sum(a);
a=a./s2;
if iter > itmax
a = rand(n);
a = a - diag(diag(a));
s1=sum(a,2);
s2=sum(a);
end
end
end
This way, I just pass in a matrix size, and a tolerance. I also use tol, instead of eps (a generally bad idea to use eps as a variable to be set), and I put in a restart for slowly convergent matrices.
Now, let me call it a few times, and see what kind of matrices tend to result.
a = bistochastic(3,eps)
a =
0 0.50635 0.49365
0.49365 0 0.50635
0.50635 0.49365 0
>> a = bistochastic(3,eps)
a =
0 0.38887 0.61113
0.61113 0 0.38887
0.38887 0.61113 0
>> a = bistochastic(3,eps)
a =
0 0.2739 0.7261
0.7261 0 0.2739
0.2739 0.7261 0
>> a = bistochastic(3,eps)
a =
0 0.28395 0.71605
0.71605 0 0.28395
0.28395 0.71605 0
>> a = bistochastic(3,eps)
a =
0 0.81008 0.18992
0.18992 0 0.81008
0.81008 0.18992 0
Do you see anything of interest? I would argue the elements all tend to be rather close to 0.5. They are not uniform at all. What does uniform mean here? To appreciate that, we need to know what the entire family of 3x3 matrices of this form is.
As it turns out, ALL doubly stochastic matrices can be written as a convex linear combination of permutation matrices. I discuss this in my answer. But since we require zero on the diagonal, we anly allow permutations that form derangements., That is, no element can remain in the original position. For a 3x3, there are only two such matrices.
a1 = [0 1 0;0 0 1;1 0 0];
a2 = [0 0 1;1 0 0;0 1 0];
We can think of those matrices as corresponding to the derangements of the vector [1 2 3], into the two vectors [2 3 1], and [3 1 2].
Now, we can write the set of all such doubly stochastic zero diagonal matrices simply as:
bistoc3 = @(t) t*a1 + (1-t)*a2;
So, for ANY value t, bistoc3 computes a 3x3 matrix of that form. Thus, if t is uniformly distributed on the interval [0,1], then bistoc3 produces a result, with elements that are also uniformly distributed in a logical sense.
If you look at the results from bistochastic above, we can see they do indeed fall into this simply family. And we can see that bistoc3 does work:
bistoc3(.1)
ans =
0 0.1 0.9
0.9 0 0.1
0.1 0.9 0
Here, we should recognize that for uniformly distributed t, we will expect to see the elements of that matrix also uniformly distributed in the interval [0,1]. As such, we should expect to see some matrices with elements near 0, and some near 1.
a13 = zeros(1,1000);
for i = 1:1000
a = bistochastic(3,1e-6);
a13(i) = a(1,3);
end
hist(a13,100)
Interesting. Instead of a uniform distribution, the elements seem to be clustered near 0.5. Remember we expect, since this is really a one-dimensional family of matrices, that the histogram should be flat.
Is that a bad thing? Well, it need not be. I mean, what is the target? Just a "random" matrix with the desired properties. Is the matrix random? Well, yes. That the distribution seems to be not uniform, is it important? Probably not terribly so, but I cannot know that for a fact.
Bruno Luong

Bruno Luong (view profile)

on 3 Sep 2019
"Remember we expect, since this is really a one-dimensional family of matrices, that the histogram should be flat. "
Not sure I would agree with John. When we introduce sum-constraint as conditional probability, even if we start with uniform distribution, it tends to have some non-uniform distribution when we look at indivisual variable. It is not possible to ensure that.
dpb

dpb (view profile)

on 3 Sep 2019
I'll have to come back and read in depth, but "I also use tol, instead of eps (a generally bad idea to use eps as a variable to be set" caught my eye. Good catch; it just slipped by that I had done that last night when threw the above together. I'll fix up my local copy and maybe repost a clean version above later on.

John D'Errico (view profile)

on 3 Sep 2019
Edited by John D'Errico

John D'Errico (view profile)

on 3 Sep 2019

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.
a =
0 0.528057118801456 0.471942881198544
0.471942881198544 0 0.528057118801456
0.528057118801456 0.471942881198544 0
a =
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]
a1 =
0 1 0
0 0 1
1 0 0
a2 = [0 0 1;1 0 0;0 1 0]
a2 =
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],
bistoc3(.3)
ans =
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),:) = []
P3 =
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),:) = []
P4 =
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,:));
end
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.
A
A =
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
>> sum(A,1)
ans =
1 1 1 1
>> sum(A,2)
ans =
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.)

Bruno Luong (view profile)

on 3 Sep 2019
Edited by Bruno Luong

Bruno Luong (view profile)

on 3 Sep 2019

For small n (up to 7/8), use RANDFIXEDSUM function on FEX
n = 5; % size of matrix
% This must be done once if n doesn't change
p = perms(1:n);
b = all(p-(1:n)~=0,2);
J = p(b,:);
m = size(J,1);
% This must be done everytime a new random matrix is requested
s = randfixedsum(m,1,1,0,1).'; % https://fr.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum
[S,I] = ndgrid(s,1:n);
R = accumarray([I(:), J(:)], S(:));
% Check
R
sum(R,1)
sum(R,2)
Result:
R =
0 0.2182 0.1291 0.4052 0.2475
0.2148 0 0.2349 0.3084 0.2419
0.2905 0.3468 0 0.1374 0.2253
0.2852 0.1408 0.2888 0 0.2853
0.2095 0.2942 0.3473 0.1490 0
ans =
1.0000 1.0000 1.0000 1.0000 1.0000
ans =
1.0000
1.0000
1.0000
1.0000
1.0000
For large n, you could adapt my asnwer in this thread, using LINPROG rather than INTLINPROG