# Random binary matrix with given rows and columns sums

I need to generate a random binary matrix $(n, n)$ whose rows sums and columns sums are $4$.

I don't manage to find a quite efficient algorithm to do this. Have you an idea please ?

NB : The distribution should be uniform.

-
Do you want it to have any particular distribution? –  joriki May 30 '13 at 9:27
If the matrix is symmetric, it's the adjacency matrix of a regular graph (with loops). You might want to take a look at algorithms for random regular graphs. –  joriki May 30 '13 at 9:33
Yes, I want the distribution to be uniform. Unfortunately, the matrix is not necessarily symmetric. –  Arnaud May 30 '13 at 9:37
Then the matrix is the adjacency matrix of a regular directed graph (with loops). I don't know if any work has been done on generating such graphs uniformly, but you might get some ideas from the algorithms for simple graphs. –  joriki May 30 '13 at 10:10

I'm not an expert in graph theory and I tried a direct method. First, I generate a matrix with the correct number of ones in the rows; then, I fix the columns.

bin <- function(n=6,k=4){ # creates a matrix w 4 ones in every row and col
A <- matrix(0,n,n) # 0 matrix
for(i in 1:n){
pos <- sample(1:n,k) # sample k positions in i-th row
A[i,pos] <- 1 # and set them to 1
}
# sumrows is now ok

sumcol <- apply(A,2,sum)
missing <- pmax(0,k-sumcol) # vector of n component, how many ones miss for each col
excess <- pmax(0,sumcol-k) # vector of n component, how many ones miss for each col
cat(sum(missing)," ",sum(excess),"\n")
#
# main idea
# pick a random column with a missing one (r) and a random row with an excess one (s)
# swap randomly a one in s-th columns with a zero in r-th column
# till matrix is ok (i.e., no columns has less than k ones)

while(sum(missing)>0){
rr <- which(excess>0)
ss <- which(missing>0)
r <- if(length(rr)==1) rr else sample(rr,1)
s <- if(length(ss)==1) ss else sample(ss,1)
ex <- which(A[,r]==1 & A[,s]==0)
j <- if(length(ex)==1) ex else sample(ex,1)
A[j,s] <- 1;A[j,r] <-0

sumcol <- apply(A,2,sum)
missing <- pmax(0,k-sumcol)
excess <- pmax(0,sumcol-k)
}
# sumcolumns is now ok

A

}


rr (ss) are vectors with the indexes of the columns with missing (excess) ones. r and s are random draw from rr and ss. there must be at least one way to randomly swap a one with a zero in columns r and s. The if statements are due to the behavior of sample in R: sample(3,1) samples 1 element from 1:3, and does not provides 3. Hence, I check for its length.

10 matrices where n=200 and k=100 are generated in about 10 seconds on my AirBook and every matrix requires on average 110 fixes (out of 100*100=10000 entries).

I'm not yet familiar with the editor, polish my text if you wish. Hope this helps.

-
I hadn't had this idea, this seems good ! The complexity is not linear unfortunately : do you think that a linear algorithm is possible ? –  Arnaud Jun 1 '13 at 21:21
On the fly, I don't know! The generation of the rows is obviously O(n). Then we fix: I expect the number of fixes to be probabilistically small as most columns cannot have a number of ones very far from k... I'll think to ways to bound (or estimate) the expected number of fixes. –  user1783444 Jun 1 '13 at 21:32
Yes, yes, I'm stupid, I believed that the operation of finding the fixes was quadratic, but it isn't... –  Arnaud Jun 1 '13 at 21:34

Expanding on Joriki's answer, you can use the algorithm inhttp://www.columbia.edu/~mo2291/Publications/DirectedRandomGraphs_29Jan13.pdf on page 7 or 8 to create a directed graph with outdegree 4 and I degree 4 at each vertex, with possible self-loops and multiple edges. Then you can exclude those with multiple edges. The directed adjacency matrix (a one in the i,j position means there is an edge from vertex i to j) gives you what you need.

The algorithm basically boils down to adding four random permutation matrices and changing the matrices if you get a 2 or higher in any entry. This is more efficient the larger $n$ is.

-
No, it doesn't work, because we cannot exclude those with multiple edges without consequences. If you are not convinced, try this configuration : 3 vertices, with outdegree 2 and indegree 2. Then link the two first vertices together, link the first vertice to itself, link the second vertice to itself. You are then forced to link the third vertice twice to itself. When you remove one of the two edges, you have the matrix $110-110-001$, which does not satisfy the conditions. –  Arnaud Jun 1 '13 at 15:19
You don't remove the multiedges, you restart the algorithm, like the monte carlo method for statistical mechanics. For small n, you may have to run it many times to get an admissible result. –  Brian Rushton Jun 1 '13 at 15:29
Not necessarily; permutation matrices are very sparse. For $n=200$, it would be highly unlikely for two permutation matrices to have a common entry. –  Brian Rushton Jun 1 '13 at 15:57
My previous deleted comment (sorry) was : for big n also, many many times ! –  Arnaud Jun 1 '13 at 15:59
I've tried to implement your method, but it always stops when we are near to the last lines. Even with $n = 2500$. –  Arnaud Jun 1 '13 at 16:01