How to solve a system of linear equations over non-negative integers?

Question

For a linear system Ax = b , where the matrix A and vector b have integer values, I want to find all non-negative integer vectors x that solve this equation.

So far, I have found some methods, such as the normal Smith form or the Erritic normal matrix form for finding integer solutions, and I assume that then I can use the linear solver to search for non-negative solutions. Is there a library that can make this easier?

Python solutions would be ideal, but if the library exists in another language, I want to know about it.

Answer 1

You can do this using integer programming by defining a non-negative integer solution variable for each available x value. Then you could solve the problem with Ax = b constraints and 0 objects that are looking for any possible integer solution for your system of equations.

This can be easily implemented using the pulp package in python. For example, consider a solution to the following system:

 x+2y+z = 5 3x+yz = 5 

You can solve this in the pulp with:

 import pulp A = [[1, 2, 1], [3, 1, -1]] b = [5, 5] mod = pulp.LpProblem('prob') vars = pulp.LpVariable.dicts('x', range(len(A[0])), lowBound=0, cat='Integer') for row, rhs in zip(A, b): mod += sum([row[i]*vars[i] for i in range(len(row))]) == rhs mod.solve() [vars[i].value() for i in range(len(A[0]))] # [1.0, 2.0, 0.0] 

This determines the feasible integer solution, x = 1, y = 2, z = 0.

Answer 2

If there exists at least one non-negative integral solution Ax=b , then integer programming should be able to find it. The set of all non-negative integral solutions can be found through the zero space A

---

Example

Using A and b in Erwin's answer:

 >>> from sympy import * >>> A = Matrix([[ 1, 2, 1], [ 3, 1,-1]]) >>> b = Matrix([20,12]) 

calculate empty space:

 >>> A.nullspace() [Matrix([ [ 3/5], [-4/5], [ 1]])] 

The empty space is one-dimensional and its basis vector is rational. Using the specific solution (2,8,2) found by Erwin, the set of all real solutions is a string parameterized as (2,8,2) + t (3,-4,5) . Since we are only interested in non-negative integral solutions, we intersect the straight line with a non-negative octant, and then with an integer lattice, which gives t ∈ {0,1,2} . Therefore, if:

• t=1 , we obtain the solution (5,4,7) .
• t=2 , we obtain the solution (8,0,12) .

Please note that these are 3 solutions found by Erwin using Z3.

---

However, it is much more complicated if the size of the zero space is greater than 1. Take a look at:

• Jesús De Loera, Many aspects of counting lattice points in polyhedra , 2005

Answer 3

Here is the Z3 model. Z3 is a proof of the theorem from Microsoft. The model is very similar to the previously proposed MIP model.

Enumerating integer solutions in MIP is not completely trivial (this can be done with some effort [link] or with the help of the “technology in advanced MIP-solvers” solution pool). With Z3 it's a little easier. It would be even simpler to use the Constraint Programming (CP) solver: they can automatically list the solutions.

Here we go:

 from z3 import * A = [[1, 2, 1], [3, 1, -1]] b = [20, 12] n = len(A[0]) # number of variables m = len(b) # number of constraints X = [ Int('x%d' % i) for i in range(n) ] s = Solver() s.add(And([ X[i] >= 0 for i in range(n) ])) for i in range(m): s.add( Sum([ A[i][j]*X[j] for j in range(n) ]) == b[i]) while s.check() == sat: print(s.model()) sol = [s.model().evaluate(X[i]) for i in range(n)] forbid = Or([X[i] != sol[i] for i in range(n)]) s.add(forbid) 

He decides

 x0+2x1+x2 = 20 3x0+x1-x2 = 12 

The solutions look like this:

 [x2 = 2, x0 = 2, x1 = 8] [x2 = 7, x1 = 4, x0 = 5] [x2 = 12, x1 = 0, x0 = 8] 

We can print the final model so that we can see how the “deny” restrictions were added:

 [And(x0 >= 0, x1 >= 0, x2 >= 0), 1*x0 + 2*x1 + 1*x2 == 20, 3*x0 + 1*x1 + -1*x2 == 12, Or(x0 != 2, x1 != 8, x2 != 2), Or(x0 != 5, x1 != 4, x2 != 7), Or(x0 != 8, x1 != 0, x2 != 12)]