Hello! As part of my efforts to develop a public toolbox for doing ecology with the WL, I am looking to optimize an often-used matrix randomization algorithm. Simply put, the end goal is an effectively random binary matrix that preserves the row and column totals of a starting binary matrix. The standard algorithm uses 'quad flips'. A 'quad' is set of four locations in a rectangular arrangement (i.e., the intersections of two rows and two columns). If the elements of the quad are either {{1,0},{0,1}} or {{0,1},{1,0}}, then the elements can be flipped while preserving the row and column totals. Randomizing a matrix thus consists of searching for such quads, and flipping them, many times. Figuring out the number of flips required to asymptotically decouple the final matrix from the starting matrix is complicated, but it rises with the size and density of the matrix and can be in the thousands or tens of thousands for some kinds matrices in ecology. And, as these randomized matrices are typically used to generate a null destruction of some statistic of interest (say, nestedness), the whole process itself may nee to be repeated maybe 1000 times. So, an efficient algorithm is key!

I attach below my first attempt, which does a single flip. For each flip it randomly picks a '1', then picks the other '1' entries in a random sequence and checks whether each forms a flippable quad (the other two entries in the quad are 0). If they do, it makes the flip and returns the new matrix, otherwise it keeps looking until every other 1 has been tested.

matrixSwap[m_] := 
 Module[{mCoords, swapCandidate, quadCandidates, n, quadCandidate},
  total = Total[m, 2];
  mCoords = Position[m, 1];
  swapCandidate = RandomSample[mCoords, 1][[1]];
  quadCandidates = RandomSample[DeleteCases[mCoords, swapCandidate]];
  n = 1;
  Until[(quadCandidate = {{m[[swapCandidate[[1]], 
         swapCandidate[[2]]]], 
       m[[swapCandidate[[1]], 
         quadCandidates[[n, 2]]]]}, {m[[quadCandidates[[n, 1]], 
         swapCandidate[[2]]]], 
       m[[quadCandidates[[n, 1]], 
         quadCandidates[[n, 2]]]]}}; (swapCandidate[[1]] =!= 
        quadCandidates[[n, 1]] && 
       swapCandidate[[2]] =!= quadCandidates[[n, 2]] && 
       Total[quadCandidate] === {1, 1} && 
       Total[quadCandidate, {2}] === {1, 1}) || n == (total - 1)), n++];
  ReplacePart[m, 
   MapThread[#1 -> #2 &, {{{swapCandidate[[1]], 
       swapCandidate[[2]]}, {swapCandidate[[1]], 
       quadCandidates[[n, 2]]}, {quadCandidates[[n, 1]], 
       swapCandidate[[2]]}, {quadCandidates[[n, 1]], 
       quadCandidates[[n, 2]]}}, Flatten[Reverse[quadCandidate]]}]]
  ]

Here is a test, using Nest to apply it many times.

testMatrix =
  RandomBinaryMatrix[100, 100, 5000];
Timing[randomizedMatrix = Nest[matrixSwap, testMatrix, 1000];]
Image[testMatrix, ImageSize -> 150]
Image[randomizedMatrix, ImageSize -> 150]
{Total[testMatrix] === Total[randomizedMatrix], 
 Total[testMatrix, {2}] === Total[randomizedMatrix, {2}]}

For me it takes about 0.8 seconds to apply 1000* swaps. Some notes on this particular algorithm:

1. My code implements a suggestion I got from Daniel Lichtblau when I asked a related question on stackexchange over 10 years ago! Basically to not generate random quads but start with pairs of ones, and to search iteratively (here using Until) for a flippable quad rather than generating all possible pairs. Good advice which I had forgotten about, but this time managed to come up more or less on my own. But, there is surely still room for improvement.
2. If more than half the entries in the matrix are 1, then it is more efficient to switch to testing pairs of zero entries. I will make this simple improvement depending on what the final version looks like.
3. The function only picks one starting '1' entry. It is possible, especially for a small matrix, that there are no flippable quads for a given starting point. That means that in n* iterations using Nest, there may be less than n actual swaps. This could be solved within the function by, if no flippable quad is found, trying again with a different starting entry. But the problem is likely not a big deal though, because small matrices are fast to randomize, so one could just make n quite large in the Nest. And large matrices are slower but almost always have quads.
4. As this is a single flip function, there is some overhead in passing the matrix, and doing some of the initial processing, e.g., the Position command, which might be avoided by writing a function that does n flips internally,

Anyway, if anyone can see some ways to significantly speed this up, I would be most interested to hear about it. A challenge perhaps? Meanwhile I will continue to see what I can do. Perhaps there is a compiled version...

Gareth