Let A be a list of n lists of m non-negative integers, such that for all j there is i with A[i][j] nonzero. Let V be a list of m positive integers.
Question: What is the fastest way to find all the lists X of n non-negative integers such that for all i then sum_j A[i][j] X[j] = V[i]?
The assumptions implies that the number of solutions is finite. See below an example of A and V, with 5499 solutions X.
Let me reformulate the problem using matrix and vector. Let A be a n-by-m matrix with non-negative integral entries and without zero column. Let V be a vector with positive integral entries. What is the fastest way to find all the vectors X, with non-negative integral entries, such that AX=V?
The usual functions for solving such a system may underuse the non-negativity. To prove so, I wrote a brute-force code finding all the solutions of such a system and applied it to an example (see below, and these crossposts on mathoverflow and on ask.sagemath), but I'm still looking for something significantly faster than this; in fact I'm looking for the fastest way.
Example
Here is the kind of system I need to solve (where x_i is non-negative integral), but with possibly more equations and variables.
[
5*x0 + 5*x1 + 5*x2 + 6*x3 + 7*x4 + 7*x5 == 24,
5*x1 + 7*x10 + 5*x6 + 5*x7 + 6*x8 + 7*x9 == 25,
5*x11 + 6*x12 + 7*x13 + 7*x14 + 5*x2 + 5*x7 == 25,
5*x12 + 6*x15 + 7*x16 + 7*x17 + 5*x3 + 5*x8 == 30,
5*x13 + 6*x16 + 7*x18 + 7*x19 + 5*x4 + 5*x9 == 35,
5*x10 + 5*x14 + 6*x17 + 7*x19 + 7*x20 + 5*x5 == 35,
5*x1 + 7*x10 + 5*x6 + 5*x7 + 6*x8 + 7*x9 == 25,
5*x21 + 5*x22 + 6*x23 + 7*x24 + 7*x25 + 5*x6 == 24,
5*x22 + 5*x26 + 6*x27 + 7*x28 + 7*x29 + 5*x7 == 25,
5*x23 + 5*x27 + 6*x30 + 7*x31 + 7*x32 + 5*x8 == 30,
5*x24 + 5*x28 + 6*x31 + 7*x33 + 7*x34 + 5*x9 == 35,
5*x10 + 5*x25 + 5*x29 + 6*x32 + 7*x34 + 7*x35 == 35,
5*x11 + 6*x12 + 7*x13 + 7*x14 + 5*x2 + 5*x7 == 25,
5*x22 + 5*x26 + 6*x27 + 7*x28 + 7*x29 + 5*x7 == 25,
5*x11 + 5*x26 + 5*x36 + 6*x37 + 7*x38 + 7*x39 == 24,
5*x12 + 5*x27 + 5*x37 + 6*x40 + 7*x41 + 7*x42 == 30,
5*x13 + 5*x28 + 5*x38 + 6*x41 + 7*x43 + 7*x44 == 35,
5*x14 + 5*x29 + 5*x39 + 6*x42 + 7*x44 + 7*x45 == 35,
5*x12 + 6*x15 + 7*x16 + 7*x17 + 5*x3 + 5*x8 == 30,
5*x23 + 5*x27 + 6*x30 + 7*x31 + 7*x32 + 5*x8 == 30,
5*x12 + 5*x27 + 5*x37 + 6*x40 + 7*x41 + 7*x42 == 30,
5*x15 + 5*x30 + 5*x40 + 6*x46 + 7*x47 + 7*x48 == 35,
5*x16 + 5*x31 + 5*x41 + 6*x47 + 7*x49 + 7*x50 == 42,
5*x17 + 5*x32 + 5*x42 + 6*x48 + 7*x50 + 7*x51 == 42,
5*x13 + 6*x16 + 7*x18 + 7*x19 + 5*x4 + 5*x9 == 35,
5*x24 + 5*x28 + 6*x31 + 7*x33 + 7*x34 + 5*x9 == 35,
5*x13 + 5*x28 + 5*x38 + 6*x41 + 7*x43 + 7*x44 == 35,
5*x16 + 5*x31 + 5*x41 + 6*x47 + 7*x49 + 7*x50 == 42,
5*x18 + 5*x33 + 5*x43 + 6*x49 + 7*x52 + 7*x53 == 48,
5*x19 + 5*x34 + 5*x44 + 6*x50 + 7*x53 + 7*x54 == 49,
5*x10 + 5*x14 + 6*x17 + 7*x19 + 7*x20 + 5*x5 == 35,
5*x10 + 5*x25 + 5*x29 + 6*x32 + 7*x34 + 7*x35 == 35,
5*x14 + 5*x29 + 5*x39 + 6*x42 + 7*x44 + 7*x45 == 35,
5*x17 + 5*x32 + 5*x42 + 6*x48 + 7*x50 + 7*x51 == 42,
5*x19 + 5*x34 + 5*x44 + 6*x50 + 7*x53 + 7*x54 == 49,
5*x20 + 5*x35 + 5*x45 + 6*x51 + 7*x54 + 7*x55 == 48
]
Here are explicit A and V from above system (in list form):
A=[
[5,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,5,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,5,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,5,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,5,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0],
[0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,7,7,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0],
[0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0],
[0,0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7]
]
V=[24,25,25,30,35,35,25,24,25,30,35,35,25,25,24,30,35,35,30,30,30,35,42,42,35,35,35,42,48,49,35,35,35,42,49,48]
Computation
I wrote a brute-force code finding all the solutions of such a system, then applied it to A, V above. It took 12 seconds to find all 5499 solutions. I'm looking for something significantly faster than this.
sage: %time LX=NonnegativeSolverPartition(A,V)
CPU times: user 11.8 s, sys: 0 ns, total: 11.8 s
Wall time: 11.8 s
sage: len(LX)
5499
Note that the time reduces to 3 seconds with PyPy3, but it is still too slow for other (bigger) such systems I need to solve.
Code
Here is my (python) code, improved by Peter Taylor (see his comment):
def NonnegativeSolverPartition(A,V):
WB = WeakBasis(A)
VB = VarBound(A,V)
PP = []
for i, ll in WB:
L = tuple(A[i][j] for j in ll)
B = tuple(VB[j] for j in ll)
PP.append(WeightedPartitionSolver(L, B, V[i]))
return list(NonnegativeSolverPartitionInter(A, V, PP, WB, [-1] * len(A[0])))
def NonnegativeSolverPartitionInter(A, V, PP, WB, X):
if any(len(P) > 1 for P in PP):
_, ii = min((len(P), i) for i, P in enumerate(PP) if len(P) > 1)
for p in PP[ii]:
PPP = PP[:ii] + [[p]] + PP[ii+1:]
Fi = Filter(PPP, list(X), WB)
if Fi:
PPPP, XX = Fi
yield from NonnegativeSolverPartitionInter(A, V, PPPP, WB, XX)
else:
assert -1 not in X
yield X
def WeakBasis(A):
return tuple(enumerate([j for j, tmp in enumerate(row) if tmp] for row in A))
def WeightedPartitions(ws, n):
def inner(ws, n):
if n == 0:
yield (0,) * len(ws)
elif ws:
w = ws[0]
lim = n // w
ws = ws[1:]
for i in range(lim + 1):
for tl in inner(ws, n - w * i):
yield (i,) + tl
return list(inner(ws, n))
def VarBound(A,V):
nvars = len(A[0])
# Solve the individual constraints and then intersect the solutions.
possible_values = [None] * nvars
row_solns = []
for row, v in zip(A, V):
lut = []
ws = []
var_assignments = []
for j, val in enumerate(row):
if val:
lut.append(j)
ws.append(val)
var_assignments.append(set())
for soln in WeightedPartitions(ws, v):
for i, w in enumerate(soln):
var_assignments[i].add(w)
for j, assignments in zip(lut, var_assignments):
if possible_values[j] is None:
possible_values[j] = assignments
else:
possible_values[j] &= assignments
return tuple(frozenset(x) for x in possible_values)
def WeightedPartitionSolver(L, B, n):
# the entries of L must be non-negative
# B gives valid values coming from other equations (see VarBound)
def inner(L, B, n):
if n == 0:
yield (0,) * len(L)
elif L:
w, allowed = L[0], B[0]
L, B = L[1:], B[1:]
for i in range(n // w + 1):
if i in allowed:
for tl in inner(L, B, n - w * i):
yield (i,) + tl
return list(inner(L, B, n))
def Filter(PP, X, W):
if [] in PP:
return None
while True:
for Wi, P in zip(W, PP):
F = FixedPoints(P)
for j in F:
P0j = P[0][j]
Wij = Wi[1][j]
if X[Wij] == -1:
X[Wij] = P0j
elif X[Wij] != P0j:
return None
LL=[]
for Wi, P in zip(W, PP):
LL.append([p for p in P if not any(X[idx] not in (-1, pval) for idx, pval in zip(Wi[1], p))])
if not LL[-1]:
return None
if PP == LL:
return LL, X
PP = LL
def FixedPoints(P):
# This would prefer P to be transposed
m=len(P)
n=len(P[0])
return tuple(i for i in range(n) if all(P[j][i] == P[0][i] for j in range(m)))
A simpler brute-force code by Max Alekseyev is available in this answer.