Find elements of an NxN matrix given constraints on row and column sums — Mathematica or Matlab

Question

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!

Answer 1

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

Answer 2

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