Ask Your Question
0

Fastest way to solve non-negative linear diophantine equations

asked 2022-06-02 21:13:15 +0100

Sébastien Palcoux gravatar image

updated 2022-06-07 12:23:08 +0100

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 flag offensive close merge delete

Comments

Use Polyhedra

FrédéricC gravatar imageFrédéricC ( 2022-06-02 21:35:04 +0100 )edit

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

Max Alekseyev gravatar imageMax Alekseyev ( 2022-06-03 03:14:54 +0100 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-03 08:18:48 +0100 )edit
add a comment

1 Answer

Sort by » oldest newest most voted
1

answered 2022-06-03 16:30:03 +0100

Max Alekseyev gravatar image

updated 2022-06-04 17:01:52 +0100

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/

edit flag offensive delete link more

Comments

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-04 05:46:24 +0100 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-04 09:07:45 +0100 )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.

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-04 09:11:37 +0100 )edit

How to use Gurobi solver?

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-04 14:05:57 +0100 )edit

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

Sébastien Palcoux gravatar imageSébastien Palcoux ( 2022-06-04 14:23:38 +0100 )edit
see more comments

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

subscribe to rss feed

Stats

Asked: 2022-06-02 21:13:15 +0100

Seen: 1,222 times

Last updated: Jun 07 '22

Related questions

Solving multilinear integer equations

How to get positive solution only of a quadratic equation

Linear diophantine equations

enumerate integer points in a polytope

Solving 3rd degree Diophantine equation in Sage

Solve symbolic matrix equations

Random positive definite matrix

Routines for Pell's equations

What is the SAGE command for calculating a Frobenius number

Plotting systems of linear equations with 3 variables

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59