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

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])

Answer 3

There is a method called block lanczos. This can your answer over a finite field. There are block lanczos solvers you find for this specific problem.

Answer 4

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 5

+1 to seberg, here is a counter example to illustrate that you should not map round:
(Sorry, it's matlab style, but you'll easily pythonize)

a =
     3     0
     0     3
     1     1
b = 
    2.71
   11.7
    0.5
x = a\b =
    0.5121
    3.5088
round(x) =
    1
    4
norm(a*round(x)-b) = 4.5193
norm(a*[0;4]-b) = 4.4367
norm(a*[1;3]-b) = 4.4299

Answer 6

I may be misunderstanding your problem, but I think you just need a combination of round and then astype(int)?

E.g.

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

x = np.linalg.lstsq(a,b)[0]
print x.round().astype(int)