How to solve a system of linear equations over the nonnegative integers?

Question

Given a linear system Ax = b, where matrix A and vector b have integer values, I want to find all nonnegative integer vectors x that solve this equation.

So far, I have found some techniques such as the Smith normal form or the Hermite normal form of a matrix to find integer solutions, and I guess I could then use a linear solver to find nonnegative solutions. Is there a library that could make this easier?

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

Answer 1

You could do this using integer programming, defining a non-negative integer decision variable for every x value you have. Then you could solve the problem with constraints Ax=b and objective 0, which searches for any feasible integer solution to your system of equations.

This can easily be implemented using the pulp package in python. For instance, consider solving the following system:

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

You could solve this in 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 identifies a feasible integer solution, x=1, y=2, z=0.

Answer 2

If at least one nonnegative integral solution of Ax=b exists, then integer programming should be able to find it. The set of all nonnegative integral solutions can be found via the null space of A.

---

Example

Using the A and b in Erwin's answer:

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

compute the null space:

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

The null space is 1-dimensional and its basis vector is rational-valued. Using the particular solution (2,8,2) that Erwin found, the set of all real solutions is a line parametrized as (2,8,2) + t (3,-4,5). Since we are interested only in nonnegative integral solutions, we intersect the line with the nonnegative octant and then with the integer lattice, which produces t ∈ {0,1,2}. Hence, if:

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

Note that these are the 3 solutions Erwin found using Z3.

---

However, it is much harder when the dimension of the null space is higher than 1. Take a look at:

• Jesús De Loera, The many aspects of counting lattice points in polytopes, 2005.

Answer 3

Here is a Z3 model. Z3 is a theorem prover from Microsoft. The model is very similar to the MIP model proposed earlier.

Enumerating integer solutions in a MIP is not totally trivial (it can be done with some effort [link] or using the "solution pool" technology in advanced MIP solvers). With Z3 this is a bit easier. Even easier could be to use a Constraint Programming (CP) solver: they can enumerate solutions automatically.

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)  

It solves

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

The solutions looks like:

[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 we can see how the "forbid" constraints 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)]