Solving linear system over integers with numpy

Question

I'm trying to solve an overdetermined linear system of equations with numpy. Currently, I'm doing something like this (as a simple example):

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

print np.linalg.lstsq(a,b)[0]

[ 1.  1.]

This works, but uses floats. Is there any way to solve the system over integers only? I've tried something along the lines of

print map(int, np.linalg.lstsq(a,b)[0])

[0, 1]

in order to convert the solution to an array of ints, expecting [1, 1], but clearly I'm missing something. Could anyone point me in the right direction?

Answer 1

You should use specialized integer problem solvers (note that integer problems are not even simple to solve). openopt is a package that for example should provide good wrappers for integer quadratic optimization such as you are doing. Trying to use linear algebra will simply not give you the correct solution that directly.

Your problem can be written as with a quadratic program, but it is an integer one, so use openopt or some other module for that. Since it is a very simple, unconstrained one, maybe there is some other approach. But for starters it is not the simple problem it looks like at first, and there are programs in openopt, etc. ready to solve this kind of thing efficiently.

Answer 2

You are looking at a system of linear diophantine equations. A quick Google search comes up with Systems of Linear Diophantine Equations by Felix Lazebnik. In that paper, the author considers the following question:

Given a system of linear equations Ax = b, where A = a(i,j) is an m × n matrix with integer entries, and b is an m × 1 column vector with integer components, does the system have an integer solution, i.e. an n × 1 solution vector x with integer components?

Answer 3

When you convert to int, the decimal part of the elements gets truncated, so it is rounded down.

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

x = np.linalg.lstsq(a,b)[0]

Result:

>>> x
array([ 1.,  1.])
>>> x[0]
0.99999999999999967
>>> x[1]
1.0000000000000002
>>> x.astype(int)
array([0, 1])
>>> map(int, x)
[0, 1]
>>> np.array([1.,1.]).astype(int) # works fine here
array([1, 1])