Solve system of linear integer equations in Python

Question

I am looking for a method to solve a system of linear equations in Python. In particular, I am looking for the smallest integer vector that is larger than all zeros and solves the given equation. For example, I have the following equation:

and want to solve .

In this case, the smallest integer vector that solves this equation is .

However, how can I determine this solution automatically? If I use scipy.optimize.nnls, like

A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])
nnls(A,b)

the result is (array([ 0., 0., 0.]), 0.0). Which is also correct, but not the desired solution...

Edit: I apologize for being imprecise in certain aspects. If anyone is interested in the details, the problem comes from the paper "Static Scheduling of Synchronous Data Flow Programs for Digital Signal Processing", Edward A. Lee and David G. Messerschmitt, IEEE Transactions on Computers, Vol. C-36, No. 1, pp. 24-35, January, 1987.

Theorem 2 says

For a connected SDF graph with s nodes and topology matrix A and with rank(A)=s-2, we can find a positive integer vector b != 0 such that Ab = 0 where 0 is the zero vector.

Directly after the proof of Theorem 2 they say

It may be desirable to solve for the smallest positive integer vector in the nullspace. To do this, reduce each rational entry in u' so that its numerator and denominator are relatively prime. Euclid's algorithm will work for this.

Answer 1

It's a wild idea perhaps but it sounds like you're looking to build a constraint solver.

minizinc is a general purpose constraint solver. Perhaps it is possible to express your constraint in such a way that minizinc can solve it?

Then it appears there is a Python library to interface with it: https://pypi.python.org/pypi/pymzn/

Answer 2

This question is quite informal as seen in the comments. Without knowing your definition of smallest (examples: l1-norm, l2-norm), it's hard to answer your specific problem.

The general problem is linked to solving a system of diophantine equations, but those are hard to solve (in the general case) and there is not much software.

A natural approach is to use Integer-programming, which is NP-hard, but commercial and some open-source solvers are very very powerful.

There is no built-in method in numpy/scipy which solves your problem without huge modifications (like: implementing some own algorithm based on numpy/scipy). Sadly, there is no IP-solver too within numpy/scipy.

Let's assume:

• smallest is the sum of variables (most of the time this is not the right approach; l2-norm is more popular, but it's a bit harder to formulate and it needs more powerful solvers)
• you want to forbid the zero-vector
• x is nonnegative

Here is some simple IP-based implementation using pulp and numpy. I don't like pulp much, but it's easy to install (pip install pulp) on all kind of popular systems.

Code

from pulp import *
import numpy as np

EPS = 1e-3

""" Input """
A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])

""" MIP """
# Variables
x = np.empty(b.shape[0], dtype=object)
for i in range(b.shape[0]):
    x[i] = LpVariable("x" + str(i), lowBound=0, upBound=None, cat='Integer')

# Problem
prob = LpProblem("prob", LpMinimize)

# Objective
prob += np.sum(x)

# Constraints
for row in range(A.shape[0]):
    prob += np.dot(A[row], x) == b[row]

prob += np.sum(x) >= EPS  # forbid zero-vector

# Solve
status = prob.solve()
print(LpStatus[status])
print([value(x_) for x_ in x])

Output

Optimal
[1.0, 1.0, 2.0]

Answer 3

Actually, it's just basic linear algebra.

>>> A = np.array([[1,-1,0], [0,2,-1],[2,0,-1]])    

Let's calculate the eigenvalues and eigenvector for this matrix.

>>> e = np.linalg.eig(A)
>>> e
(array([ -4.14213562e-01,  -1.05471187e-15,   2.41421356e+00]), array([[-0.26120387, -0.40824829,  0.54691816],
       [-0.36939806, -0.40824829, -0.77345908],
       [-0.89180581, -0.81649658,  0.32037724]]))    
>>>> np.round(e[0], 10)
array([-0.41421356, -0.        ,  2.41421356])

After rounding, we see that 0 is an eigenvalue of our matrix A. So, the eigenvector s for the 0-eigenvalue is a good candidate for your equation system.

>>> s = e[1][:,1]
>>> s
array([-0.40824829, -0.40824829, -0.81649658])    

Multiplying an eigenvector with the related matrix results in the eigenvector itself, scaled by the related eigenvalue. So, in the above case we see: As = 0s = 0

>>> np.round(A.dot(s), 10)
array([ 0.,  0.,  0.])

Since we are interested in an integer solution, we have to scale the solution vector.

>>> x = s / s[1]
>>> x
array([ 1.,  1.,  2.])

Hope this answer solves your problem.

Answer 4

To find the exact solution that you want, numpy and scipy are probably not the best tools. Their algorithms generally work in floating point, and are not guaranteed to give the exact answer.

You can use sympy to get the exact answer to this problem. The Matrix class in sympy provides the method nullspace() that returns a list of basis vectors for the nullspace. Here's an example:

In [20]: from sympy import Matrix, lcm

In [21]: A = Matrix([[1, -1, 0], [0, 2, -1], [2, 0, -1]])

Get the vector in the nullspace. nullspace() returns a list; this code assumes that the rank of A is 2, so the list has length one:

In [22]: v = A.nullspace()[0]

In [23]: v
Out[23]: 
Matrix([
[1/2],
[1/2],
[  1]])

Find the least common multiple of the denominators of the values in v so that we can scale the result to the smallest integers:

In [24]: m = lcm([val.q for val in v])

In [25]: m
Out[25]: 2

x holds the integer answer:

In [26]: x = m*v

In [27]: x
Out[27]: 
Matrix([
[1],
[1],
[2]])

To convert that result to a numpy array of integers, you can do something like:

In [52]: np.array([int(val) for val in x])
Out[52]: array([1, 1, 2])

Answer 5

There is a solver for this kind of equation at pypi. It apparently can compute the Hermite normal form of a matrix which in turn can be used to solve your problem.

The package is versioned 0.1.

Sage also appears to support the Hermite normal form.

The special case of a homogeneous system, i.e. b=0 is a bit easier, here is a solver for the simplest possible case of the matrix having rank n-1

import sympy
import numpy as np

def create_rd(n, defect=1, range=10):
    while True:
        res = sympy.Matrix((np.random.randint(-range+1,range,(n, n-defect))
                           @np.random.randint(0,2,(n-defect, n)))
                           .astype(object))
        if res.rank() == n-defect:
            break
    return res

def solve(M):
    ns = M.nullspace()
    ns = [n / sympy.gcd(list(n)) for n in ns]
    nsnp = np.array([[int(k) for k in n] for n in ns])
    if len(ns) == 1:
        return ns[0], nsnp[0]
    else:
        raise NotImplementedError

Sample output:

>>> M = create_rd(4)  # creates a rank-deficient matirx
>>> ns, nn = solve(M) # finds the 1d nullspace and a minimal integer basis vector
>>> M
Matrix([
[-7, -7, -7, -12],
[ 0,  6,  0,   6],
[ 4,  1,  4,  -3],
[-4, -7, -4,  -9]])
>>> ns
Matrix([
[-1],
[ 0],
[ 1],
[ 0]])
>>> M*ns
Matrix([
[0],
[0],
[0],
[0]])
>>> M = create_rd(40) # we can go to higher dimensions
>>> ns, nn = solve(M) # but solutions quickly become unwieldy
>>> ns
Matrix([
[  8620150337],
[-48574455644],
[ -6216916999],
[-14703127270],
 < - snip - >