One possibility to consider is generating ab initio given the row sums and column sums. The code below will do this except it is not guaranteed to terminate (oh well). The idea is to generate rows that conform to the desired sums, then swap 0-1 in column pairs that have a deficit and an excess of 1s respectively, choosing a row with such a 0-1 pair. We put in place a stop after iterating a certain number of times so we don't get painted into a corner so to speak (this appens if for some pair of columns as described, there is no row with the requisite 0-1 pair to swap). When this happens we restart the process.
genBinaryMatrix[rsums_, csums_] := Catch[Module[
{nr = Length[rsums], nc = Length[csums], mat, posns, cdiffs,
locols, hicols, done, clow, chigh, rrow},
While[True,
mat = Normal[SparseArray[Flatten@Table[
posns = RandomSample[Range[nc], rsums[[i]]];
Map[{i, #} -> 1 &, posns], {i, nr}]]];
cdiffs = Total[mat] - csums;
locols = Flatten@Position[cdiffs, _?Negative];
hicols = Flatten@Position[cdiffs, _?Positive];
If[(locols == {} && hicols == {}), Throw[mat, "done"]];
done = (locols == {} && hicols == {});
jj = 0;
While[jj < nr*nc,
jj++;
clow = RandomChoice[locols];
chigh = RandomChoice[hicols];
rrow = RandomChoice[Flatten[Position[mat[[All, clow]], 0]]];
If[mat[[rrow, chigh]] == 0, Continue[]];
mat[[rrow, clow]] = 1;
mat[[rrow, chigh]] = 0;
cdiffs[[clow]]++;
If[cdiffs[[clow]] == 0,
locols = DeleteElements[locols, {clow}]];
cdiffs[[chigh]]--;
If[cdiffs[[chigh]] == 0,
hicols = DeleteElements[hicols, {chigh}]];
done = (locols == {} && hicols == {});
If[(locols == {} && hicols == {}), Throw[mat, "done"]];
];
]], "done"]
Here is an experiment. We create a matrix, check the row and column sums, and generate a new one based thereon.
SeedRandom[1234];
mat = RandomInteger[1, {8, 8}]
rowsums = Map[Total, mat]
colsums = Total[mat]
newmat = genBinaryMatrix[rowsums, colsums]
Map[Total, newmat] == rowsums
Total[newmat] == colsums
(* Out[172]= {{1, 1, 1, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 1, 0, 0, 1}, {0, 1,
1, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 0, 1, 1, 1}, {0, 0, 0, 0, 0, 0, 0,
0}, {1, 0, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 0, 1, 0, 0, 0}, {0, 1, 0, 1,
1, 0, 0, 1}}
Out[173]= {3, 4, 4, 3, 0, 2, 2, 4}
Out[174]= {2, 5, 3, 1, 4, 2, 2, 3}
Out[175]= {{0, 0, 1, 0, 0, 1, 0, 1}, {0, 1, 1, 0, 1, 0, 0, 1}, {1, 1,
0, 1, 0, 0, 1, 0}, {1, 1, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0,
0}, {0, 1, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 0, 0}, {0, 1, 0, 0,
1, 0, 1, 1}}
Out[176]= True
Out[177]= True *)
To give an idea of speed, I show this for a 200 x 200.
SeedRandom[1234];
dims = 200;
mat = RandomInteger[1, {dims, dims}];
rowsums = Map[Total, mat];
colsums = Total[mat];
Timing[newmat = genBinaryMatrix[rowsums, colsums];]
Map[Total, newmat] == rowsums
Total[newmat] == colsums
(* Out[183]= {1.42914, Null}
Out[184]= True
Out[185]= True *)