2

I'm trying to find the elements of an NxN matrix so that given row and column sums are satisfied. The only other constraint is that the elements of the matrix are restricted to the integers 0 and 1.

This is probably a few simple lines in Matlab or Mathematica (I've hunted around on here and found a few related questions in R), but my relevant coding experience is only in constrained optimization -- I'm lost here without an objective function to minimize.

My first thought was to set it up a linear algebra problem, Ax = b, but that doesn't work because A'A is singular. (I'll provide detail on this in a comment below just in case it helps.)

My next thought was to try to "trick" the NMinimize solver in Mathematica by giving it a "dummy" objective function subject to the constraint Ax = b and the integer constraints:

    x = Map[Subscript[x, #] &, Range[1, Ngrid^2]];
    x = Map[Subscript[x, #] &, Range[1, Ngrid^2]];
    NMinimize[{x.x\[Transpose], A.x == b && Apply[And,Thread[GreaterEqual[x, 0]]] && 
    Apply[And, Thread[LessEqual[x, 1]]] && Element[x,Integers]}, x, 
    Method -> "DifferentialEvolution", MaxIterations -> 2000];

But it encounters the same "recursion" error.

I have a feeling this isn't the way to approach the problem. Is there some command that allows me to populate elements of a matrix according to constraints? If it doesn't generate a unique solution, is there a way to view the set of solutions?

Thanks, Dan

EDIT: Here is what I had in mind for the linear algebra solution, Ax = b:

A is a 2N x N^2 rectangular matrix of 0s and 1s. Each column corresponds to one of the N^2 elements in the solution vector x (the reshaped NxN matrix we're after). (EDIT: The first N rows correspond to the elements of x that are affected by each of the row constraints, the second N rows correspond to the elements of x that are affected by each of the column constraints.)

x is an N^2 x 1 vector of 0s and 1s (the reshaped NxN solution matrix).

b is a 2N x 1 vector of the row and column sums (the constraints).

Example:

If N=3 and we want to find x (ETA: The "x" below is the solution we're after; we don't know it ahead of time):

    1   0   0
    1   0   1
    0   1   0

We could write A = 

    1   1   1   0   0   0   0   0   0
    0   0   0   1   1   1   0   0   0
    0   0   0   0   0   0   1   1   1
    1   0   0   1   0   0   1   0   0
    0   1   0   0   1   0   0   1   0
    0   0   1   0   0   1   0   0   1

And z = 

    1
    2
    1
    2
    1
    1

Solving isn't as simple as x = inv(A'A)*A'z, because A'A doesn't have an inverse.

Thanks in advance for any help you can provide!

DanMcK
  • 37
  • 9
  • 1
    Are you trying to find the `sum` of rows and columns that match a certain value? Or are you trying to find a subset of the `A` matrix that would create a new `NxN` matrix with the proper values? Really confused by what you want... – Matt Dec 18 '15 at 15:30
  • Also, should the sums be different for each row and column? – Dr. belisarius Dec 18 '15 at 15:39
  • Hi all, sorry for the confusion. I'm trying to find an NxN matrix so that constraints on row and column sums are satisfied. The sums can be different for each row/column -- these are given, and I have to find the NxN matrix (with the restriction that each element of the matrix is 0 or 1). I think the confusion is coming from my example -- assume you're only given "b" (which tells you the first row must sum to 1, the second row to 2, ... , the first column to 2, and so on) and you know x must be 3x3, using only integers 0 and 1. How could you find x? Thanks again. – DanMcK Dec 18 '15 at 16:26

2 Answers2

6

This is a naive approach. It depends on how you need it to scale up:

n = 3;
sumsrows = {2, 1, 1};
sumscols = {1, 2, 1};
X = Array[x, {n, n}];

fi = FindInstance[
   And @@ Thread[(Tr /@ X) == sumsrows] && 
   And @@ Thread[(Tr /@ Transpose@X) == sumscols] && 
   And @@ Thread[0 <= Flatten@X <= 1], Flatten@X, Integers]

Partition[fi[[1, All, 2]], n] // MatrixForm

Mathematica graphics

Dr. belisarius
  • 60,527
  • 15
  • 115
  • 190
  • Thank you very much @Dr. belisarius -- this is exactly what I was looking for. I'll tinker with some larger N and see how it holds up. I really appreciate the help. – DanMcK Dec 18 '15 at 18:02
  • Is there a way to adapt this to the case of a rectangular matrix? – Semiclassical Aug 04 '17 at 16:53
0

Another very naive approach would be to generate random "binary" matrices and stop if one of them meets your criteria. Here is the MATLAB code:

n = 3;
row_sums = [2; 1; 1];
col_sums = [1, 2, 1];
solution_found = false;

while ~solution_found
    X = randi([0, 1], n);
    if isequal(sum(X, 2), row_sums) && isequal(sum(X, 1), col_sums)
        solution_found = true;
    end
end
statmerkur
  • 429
  • 3
  • 16