3
|
471 Views
|
7 Replies
|
8 Total Likes
View groups...
Share
GROUPS:

Optimization challenge: row and column total preserving matrix randomization

Posted 19 days ago
7 Replies
Sort By:
Posted 12 days ago
 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 *)
Posted 12 days ago
 Interesting! By starting with random elements within rows, and then adjusting them to create the necessary column sums, the resultant matrix will presumably be fully random within the row and column sum constraints. Therefore one can stop as soon as the constraints are achieved, rather than 'swapping' for a pre-determined number of times to move sufficiently far away from the original matrix, as in the algorithm I am using.In practice, where my code takes 1.2 seconds to do 10,000 swaps, yours takes an average of 0.8 seconds to achieve the constraints, but with an exponential distribution of actual times. (The longest time was a bit over 4 seconds.)I need to check if the long runs might have failed to meet the constraints... ;)
Posted 13 days ago
 What are the requirements of the result? Need it be a succession of these quad flips from a given matrix? Or is it acceptable simply to produce a new binary matrix with prescribed row and column sums?
Posted 13 days ago
 Just a new matrix. I'm attaching a notebook with what I have so far. It is the MatrixSwap function that is at play here: the MatrixRandomize function calls that (many times) when the option to preserve row and column totals is selected. (Otherwise it calls various simpler functions which don't really need to be optimized further.) I implement the refinement to find quads based on the less common element (1 or 0).
Posted 15 days ago
 You might get a big speed boost if you use Compile and just rewrite in a procedural programming style.As for the advice you received, I hope it was good. I’m not sure I’d have thought of it today (and certainly didn’t remember it).
Posted 15 days ago
 I’m still thinking about your code and I’m not ready to suggest a new algorithm. But in any case speed will increase if instead of ordinary arrays use SparseArrays.
Posted 13 days ago
 I think you are right in a sense: note that in the MatrixSwap function I pull the coordinates of the target elements (1 or 0) using Position, and thereafter work with that list, which is effectively a Sparse Array. (Although I look back at the original matrix for the 'opposite corners'. Maybe it would be faster to do that with a Sparse Array, but I suspect not.)