Generate random Bistochastic matrix with zero diagonals

3 次查看(过去 30 天)
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.

采纳的回答

dpb
dpb 2019-9-3
编辑:dpb 2019-9-3
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.
  11 个评论
Bruno Luong
Bruno Luong 2019-9-3
编辑:Bruno Luong 2019-9-3
"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 2019-9-3
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.

请先登录,再进行评论。

更多回答(2 个)

John D'Errico
John D'Errico 2019-9-3
编辑:John D'Errico 2019-9-3
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
Bruno Luong 2019-9-3
编辑:Bruno Luong 2019-9-3
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

类别

Help CenterFile Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

产品


版本

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by