# Finding all possible $n\times n$ matrices with non-negative entries and given row and column sums.

What I am ultimately attempting to do is find the solution that maximizes a given equation so I need to find all possible solutions so I can check them. I need all possible solutions to an $n \times n$ matrix, given the row and column totals. All values must be non-negative. Other than that the row and sum columns are the only constraints.

• Are you asking about (non-negative) integer entries? – hardmath Feb 27 '13 at 16:00
• If you are, this is the problem of generating all feasible contingency tables for the specified row and column sums. – hardmath Feb 27 '13 at 16:06
• Yes, they are non-negative integers. Sorry, I should have specified. This seems like the problem I am encountering. – Simon Feb 27 '13 at 16:23
• Take a look at the answer I gave on this previous Question by Douglas S. Stones. Although the framing is in terms of multisets, we sorted out that it is mathematically equivalent. So, two notes: there are typically a large number of these, and counting them is a "known hard problem" (in the sense of #P-completeness). It would probably be constructive to focus on what kind of objective function you are working with, to see if it allows some shortcuts in the search space (symmetries, etc.). – hardmath Feb 27 '13 at 16:29
• See also the link given by @TmothyWagner for a simpler Question to a paper Matrices with Prescribed Row and Column Sums by Alexander Barvinok (2010), a couple of years more recent than the survey of the latter author in my answer. – hardmath Feb 27 '13 at 16:53

It might be interesting to mention the RSK correspondence (although it might not be of sufficient help if $n$ is really very large). It gives an algorithmic bijection between the set of non-negative integer matrices with column and rows sums $\alpha,\beta$ respectively (these are vectors of non-negative integers) and pairs $(P,Q)$ of equal shape semistandard Young tableaux of respective weights $\alpha,\beta$. Therefore the number of matrices is given by $$\sum_\lambda K_{\lambda,\alpha}K_{\lambda,\beta}$$ where $\lambda$ is a partition of $\sum\alpha=\sum\beta$, and $K_{\lambda,\alpha}$ is the number of semistandard Young tableaux of shape $\lambda$ and weight $\alpha$, which is a Kostka number. This at least gives a method to count the number of matrices that is much more efficient than enumerating them all (because of the product in the formula and the fact that the $K_{\lambda,\alpha}$ are non-negative). There are no really simple formulas for computing Kostka numbers, but since they are weight multiplicities, there is a recursive formula that computes them again much more efficiently than by counting semistandard Young tableaux.

In principle there is an "efficient" recursive procedure for generating all the possible non-negative integer $n \times n$ matrices with prescribed row and column sums (contingency tables), but we should be forewarned that the number of such solutions may be both enormous and difficult to estimate with precision.

The simpler (but still voluminous) generation of $2 \times n$ contingency tables is the subject of this earlier Question, and a more general problem with (essentially) $m \times k$ matrices was discussed here, including an informally stated recursive generating procedure.

Often authors contemplate the problem of randomly selecting solutions with uniform probability, but here the stated goal is an exhaustive search/brute-force determination of a maximized objective function. If more were known about the structure of that objective function, likely some practical shortcuts could be found (e.g. exploiting symmetries of the search space).

The Question asks us to generate all possible matrices with prescribed row and column sums. Counting them need not concern us directly. Any algorithm that avoids duplication or unnecessary trial-and-error is “efficient”. Below we present Prolog code implementing such an algorithm, predicate contingencyTable/3, taking as input two lists, of row sums and column sums respectively, and producing on backtracking all possible non-negative integer matrices satisfying those constraints.

balanceCondition/2 Check that totals of row sums and column sums agree, a necessary and sufficient condition for solutions to exist.

generateTable/3 Select non-negative entries for the first row (see generateRow) to match its row sum and reduce the respective column sums (but not to less than zero). The “balance condition” is preserved, so the procedure can recursively fill the next row. When one row is left, selection there is unique, and at that point the sum of the last reduced column sums equals the last row sum.

generateRow/4 Entries for one row can also be selected recursively. The first entry in a row will have two lower bounds and two upper bounds. While one lower bound is trivial (non-negativity), the other is more subtle: assign enough that later column sums do not limit the row sum from being attained. Both upper bounds are obvious; the entry cannot exceed the row sum or the column sum when the entry is selected. Once an entry within the allowable bounds is chosen, reduce the row sum accordingly and apply the same procedure to the rest of the row. The last row entry is uniquely determined.

/*
Generate 2-way contingency tables with specified margins,
non-negative matrices with prescribed row & column sums.
*/

/* contingencyTable(+RowSums,+ColSums,?Table) */
contingencyTable(RowSums,ColSums,Table) :-
balanceCondition(RowSums,ColSums),       /* feasible? */
generateTable(RowSums,ColSums,Table).

balanceCondition(RowSums,ColSums) :-
sum(RowSums,0,Total),
sum(ColSums,0,Total).

generateTable([_],ColSums,[ColSums]) :- !.
generateTable([RowSum|RowSums],ColSums,[Row|Rows]) :-
generateRow(RowSum,ColSums,ReducedColSums,Row),
generateTable(RowSums,ReducedColSums,Rows).

generateRow(RowSum,[ColSum],[ColRed],[RowSum]) :-
!,
ColRed is ColSum - RowSum.
generateRow(RowSum,[ColSum|ColSums],[ColRed|ColReds],[R|Rs]) :-
sum(ColSums,0,Sum),
(Sum > RowSum -> RMin = 0 ; RMin is RowSum - Sum),
(ColSum > RowSum -> RMax = RowSum ; RMax = ColSum),
between(RMin,RMax,R),
RowRed is RowSum - R,
ColRed is ColSum - R,
generateRow(RowRed,ColSums,ColReds,Rs).

/* helper predicates sum/3 and between/3 */
sum([ ],Sum,Sum).
sum([H|T],Start,Total) :-
SoFar is Start + H,
sum(T,SoFar,Total).

/* argument order compatible with SWI-Prolog builtin */
between(X,Z,X) :- X =< Z.
between(X,Z,_) :- X > Z, !, fail.
between(X,Z,K) :-
Y is X+1,
between(Y,Z,K).