# Fastest way to solve non-negative linear diophantine equations

Let $A$ be a matrix in $M_{n \times m}(\mathbb{Z}_{\ge 0})$ without zero column.

Let $V$ be a vector in $\mathbb{Z}_{> 0}^m$.

What is the fastest way to find all the solutions $X \in \mathbb{Z}_{\ge 0}^n$ of the following equation? $$AX=V$$ By non-negativity and "without zero column" assumption, there always are finitely many solutions only.

The usual functions on SageMath for solving such a system may underuse the non-negativity. To prove so, I wrote a naive code finding all the solutions of such a system (see below), but I'm still looking for something significantly faster than it (see this crosspost on mathoverflow).

Here is the kind of system I need to solve (where $x_i \in \mathbb{Z}_{\ge 0}$), 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]

I wrote a naive code finding all the solutions of such a system, then applied it to A, V above. It took 41 seconds to find all 5499 solutions. I'm looking for something significantly faster than this.

sage: %time LX=NonnegativeSolverPartition(A,V)
CPU times: user 40.7 s, sys: 0 ns, total: 40.7 s
Wall time: 40.9 s
sage: len(LX)
5499


Here is my code:

import copy

def NonnegativeSolverPartition(A,V):
global LX
WB=WeakBasis(A)
VB=VarBound(A,V)
m=len(A[0])
PP=[]
LX=[]
Or=[]
c=0
for l in WB:
[i,ll]=l
L=[A[i][j] for j in ll]
B=[VB[j] for j in ll]
n=V[i]
P=WeightedPartitionSolver(L,B,n)
Or.append([len(P),c])
c+=1
PP.append(P)
Or.sort()
for i in range(len(Or)):
if Or[i][0]!=1:
break
ii=Or[i][1]
for p in PP[ii]:
X=[-1 for i in range(m)]
PPP=copy.deepcopy(PP)
PPP[ii]=[p]
Fi=Filter(PPP,X,WB)
if Fi!=[]:
[PPPP,XX]=Fi
NonnegativeSolverPartitionInter(A,V,PPPP,WB,VB,XX)
return LX

def NonnegativeSolverPartitionInter(A,V,PP,WB,VB,X):
global LX
Or=[[len(PP[i]),i] for i in range(len(PP))]
Or.sort()
c=0
for i in range(len(Or)):
if Or[i][0]!=1:
c=1
break
ii=Or[i][1]
if c==0:
if -1 in X:
print('problem: -1 in X')
else:
if V==[sum([A[i][j]*X[j] for j in range(len(A[0]))]) for i in range(len(A))]:
LX.append(X)
else:
for p in PP[ii]:
PPP=copy.deepcopy(PP)
PPP[ii]=[p]
Fi=Filter(PPP,X,WB)
if Fi!=[]:
[PPPP,XX]=Fi
NonnegativeSolverPartitionInter(A,V,PPPP,WB,VB,XX)

def WeakBasis(A):
n=len(A)
m=len(A[0])
L=[]; F=[]
for i in range(n):
FF=[]
for j in range(m):
if A[i][j]!=0:
FF.append(j)
F.append(FF)
for j in range(m):
LL=[]
for i in range(n):
if A[i][j]!=0:
LL.append([len(F[i]),i])
else:
LL.append([0,i])
LL.sort(reverse=true)
L.append(LL[0])
return [[i,F[i]] for i in range(n)]

def VarBound(A,V):
M=[]
n=len(A)
m=len(A[0])
for j in range(m):
N=[]
for i in range(n):
if A[i][j]!=0:
N.append(V[i]//A[i][j])
M.append(min(N)+1)
return M

def WeightedPartitionSolver(L,B,n): # the entries of L must be non-negative
global P        #B for upperbound coming from other equations (see VarBound)
P=[]
l=len(L)
S=[0 for i in range(l)]
m=n//L[0]+1
for i in range(min(m,B[0])):
S[0]=i
WeightedPartitionSolverInter(L,B,n-L[0]*i,1,S)
return P

def WeightedPartitionSolverInter(L,B,n,j,S):
global P
l=len(L)
if j==l:
if n==0:
P.append(S)
else:
if n==0:
P.append(S)
else:
m=n//L[j]+1
for i in range(min(m,B[j])):
SS=copy.deepcopy(S)
SS[j]=i
WeightedPartitionSolverInter(L,B,n-L[j]*i,j+1,SS)

def Filter(PP,X,W):
if [] in PP:
return []
XX=copy.deepcopy(X)
LL=[]
for i in range(len(PP)):
P=PP[i]
F=FixedPoints(P)
Wi=W[i]
for j in F:
P0j=P[0][j]
Wij=Wi[1][j]
if XX[Wij]==-1:
XX[Wij]=P0j
elif XX[Wij]!=P0j:
return []
for i in range(len(PP)):
L=[]
P=PP[i]
Wi=W[i]
lP0=len(P[0])
for p in P:
c=0
for j in range(lP0):
Wij=Wi[1][j]
if XX[Wij]!=-1:
if p[j]!=XX[Wij]:
c=1
break
if c==0:
L.append(p)
if L==[]:
return []
LL.append(L)
if [PP,X]==[LL,XX]:
return  [LL,XX]
else:
return Filter(LL,XX,W)

def FixedPoints(P):
m=len(P)
n=len(P[0])
if m==1:
return [i for i in range(n)]
F=[]
for i in range(n):
c=0
for j in range(m):
if P[j][i]!=P[0][i]:
c=1
break
if c==0:
F.append(i)
return F

edit retag close merge delete

Use Polyhedra

( 2022-06-02 21:35:04 +0200 )edit

I suggest to employ normaliz backend and use .integral_points()function. See https://doc.sagemath.org/html/en/refe...

( 2022-06-03 03:14:54 +0200 )edit

@Max Alekseyev Can you show me on above example? I added the explicit matrix $A$ and vector $V$ (in list form).

( 2022-06-03 08:18:48 +0200 )edit

Sort by » oldest newest most voted

There are two essential problems to consider here:

• (i) finding ONE solution or establishing that there are NONE;
• (ii) finding ALL solutions.

In my experience, problem (i) is best solved with (M)ILP solvers, while problem (ii) is solved with polyhedral computation, which may be significantly slower than ILP. So, even in problem (ii) if there is a possibility for no solutions, it's worth to run ILP first to rule out this possibility, and only when one solution is found proceed with finding them all.

from sage.numerical.mip import MIPSolverException

def mysolve(A,V):
n = len(A)      # number of equations
assert n == len(V)
m = len(A[0])   # number of variables

milp = MixedIntegerLinearProgram(solver='ppl')
milp.set_objective( None )
x = milp.new_variable(integer=True, nonnegative=True)
for i in range(n):
milp.add_constraint( sum(A[i][j] * x[j] for j in range(m)) == V[i] )

try:
milp.solve()
except MIPSolverException:
return []       # no solutions

X = milp.get_values(x)
print(f'one solution is found: {[X[i] for i in range(m)]}')
print('proceed with finding all solutions')

# --------------------------------------------------------------

# we can extract polyhedron from MILP, or construct it explicitly
# P = milp.polyhedron( backend='normaliz', base_ring=QQ )

# explicit polyhedron construction:
equalities = [ [-V[i]] + A[i] for i in range(n) ]
inequalities = [ [0]*(j+1) + [1] + [0]*(m-j-1) for j in range(m) ]
P = Polyhedron( backend='normaliz', ieqs=inequalities, eqns=equalities, base_ring=QQ )
print('Polyhedron:', P)

return P.integral_points()


In the given example, one solution is found momentarily (tested with ppl and Gurobi solvers), but I was not patient enough to get them all.

PS. Notice that for a polyhedron extracted from MILP, it may be tricky to map the coordinates of its integral points back to the variables - see https://ask.sagemath.org/question/55520/

more

How long do your computer need to get 'one solution is found'?

( 2022-06-04 05:46:24 +0200 )edit

I did not get it on my laptop after 1 hour.

( 2022-06-04 09:07:45 +0200 )edit

I think that this approach is under exploiting the fact that the entries of $A$, $V$ and $X$ are non-negative integers. It may be the fastest way for integers but not for non-negative integers.

( 2022-06-04 09:11:37 +0200 )edit

How to use Gurobi solver?

( 2022-06-04 14:05:57 +0200 )edit

I expect that a code really exploiting the non-negative integrality should find all the solutions within a second.

( 2022-06-04 14:23:38 +0200 )edit