Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am implementing the 'swap' algorithm for a binary matrix to generate matrix permutations that maintain row and column totals. My problem is that my matrices are large and sparse (e.g., 19774 x 942, 0.14% fill). The standard algorithm is this:

Pick two rows and two columns at random, and examine the 2x2 sub matrix they define. If it is either {{1,0},{0,1}} or {{0,1},{1,0}}, then flip the entries. Otherwise do nothing and try another combination of rows and columns. Here is my version, which is not at all optimized:

 swap[matrix_] := 
 Module[{r, c, subMatrix, rowsChoice, colsChoice, outout},
  {r, c} = Dimensions[matrix];
  subMatrix = {0, 0, 0, 0};
  While[(subMatrix =!= {1, 0, 0, 1}) && (subMatrix =!= {0, 1, 1, 0}),
   rowsChoice = RandomInteger[{1, r}, 2];
   colsChoice = RandomInteger[{1, c}, 2];
   subMatrix = 
    Extract[matrix, 
     i = Flatten[Outer[List, rowsChoice, colChoice], 1]];       
   ];
  output = 
   ReplacePart[matrix, 
    MapThread[#1 -> #2 &, {i, 
      If[subMatrix == {1, 0, 0, 1}, {0, 1, 1, 0}, {1, 0, 0, 1}]}]]
  ]

Optimizations aside, obviously, for a large sparse matrix, the overwhelming majority of the random submatrices do not satisfy the swap criterion (most are {{0,0},{0,0}}), and so the process is quite slow.

So, I'm trying to come up with an efficient algorithm that does not just shoot in the dark. Some simple possibilities are to constrain the row/column choice, for example by making sure that there are at least SOME ones in the submatrix. Another idea I'm working on is to first make a list of all the submatrices that satisfy the 'swappable' condition. Then one simply chooses one of these. Making the first list is fairly easy, and not too slow given that the matrix is sparse

allSubMatrices[matrix_] := 
 Module[{ones, fill, oneCorner, otherOneCorners, oneCorners, 
   otherCorners, cornerFill},
  ones = Position[matrix, 1];
  fill = Length[ones];
  Cases[Flatten[Table[
     oneCorner = ones[[i]];
     otherOneCorners = 
      Select[Drop[ones, 
        i], #[[1]] =!= oneCorner[[1]] && #[[2]] =!= 
          oneCorner[[2]] &];
     Table[
      oneCorners = {oneCorner, otherOneCorners[[j]]};
      otherCorners = {{oneCorners[[1, 1]], 
         oneCorners[[2, 2]]}, {oneCorners[[2, 1]], 
         oneCorners[[1, 2]]}};
      cornerFill = Extract[matrix, otherCorners];
      If[cornerFill == {0, 0}, {{oneCorners[[1, 1]], 
         oneCorners[[2, 1]]}, {oneCorners[[1, 2]], 
         oneCorners[[2, 2]]}}, {}], {j, 1, 
       Length[otherOneCorners]}], {i, 1, 1}], 1], {{_, _}, {_, _}}]
  ]

The output is a list pairs of rows and columns that define 'swappable' submatrices. The harder trick is, after doing a swap, to modify this list to reflect the new set of submatrices. A swap changes four matrix entries, and removes all valid submatrices that depended on these entries. Those are relatively easy to find in the existing list. But a swap potentially adds a number of newly valid submatrices, which need to found as well, and I haven't yet figured out how to find those quickly.

Before I go further down this rabbit hole, I'm wondering if the stackexchange geniuses have any suggestions about the best way to speed up the swap algorithm. Should I stick with the basic algorithm and just work on optimizing/compile it? Is there a known (or obvious) more efficient algorithm for the sparse case that I'm missing?

Thanks,

Gareth

share|improve this question
    
Thanks for the suggestion. But all those cases deal with contiguous submatrices, whereas I am also looking for non-contiguous ones -- the intersections of any two rows and any two columns. But I'll see if there are any hints in them that I can use. –  Gareth Dec 21 '13 at 21:09
    
Since the number of 0s and 1's remains fixed you could have a table of the positions of the 1s, updated after each swap. Then you can select random pairs from there, say positions {a,b} and {c,d}. Then check whether the two associated positions {a,d} and {c,b} are 0s. If so then do a swap and update the 1s position list. Since the matrix is sparse this should "usually" lead to swaps rather than wasted attempts. –  Daniel Lichtblau Dec 21 '13 at 21:10
    
@Gareth You're right, I didn't think carefully before posting. Sorry about that. Since it's not relevant I'll just delete the comment, and later this one too. –  Szabolcs Dec 22 '13 at 0:49
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.