7

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?

arshajii
  • 127,459
  • 24
  • 238
  • 287
  • 3
    You really should clarify that question. The text states overdetermined, so maybe I am interpreting it wrong, but that suggests to me that you do not know that there is an exact solution, also I don't see a reason why the best/exact solution should be integer (if you know that one of these its a huge missing information). In which case the example oversimplifies though, because it has one best (even exact) integer solution. – seberg Dec 16 '12 at 13:39

8 Answers8

7

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.

seberg
  • 8,785
  • 2
  • 31
  • 30
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?

NPE
  • 486,780
  • 108
  • 951
  • 1,012
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])
Akavall
  • 82,592
  • 51
  • 207
  • 251
  • I figured it was something like this - so is there any way to use only integers in the computations to avoid these float inaccuracies? – arshajii Dec 16 '12 at 03:19
  • I think `np.round` should do it. – Akavall Dec 16 '12 at 03:20
  • 1
    Sorry, but this happens to work here, because the exact solution happens to be integer for this example. There is no reason rounding/truncating will give you the best integer solution. – seberg Dec 16 '12 at 11:29
1

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)
Joe Kington
  • 275,208
  • 71
  • 604
  • 463
  • 1
    Thanks, this also helps a lot. So numpy doesn't have a function specifically for solving linear systems over integers only? I was planning on using a rounding method like this as a last resort, but it shouldn't make too much of a difference. – arshajii Dec 16 '12 at 03:34
  • 3
    @A.R.S.: do you want the linear algebra to be solved using integer arithmetics or do you want the result to be of integer datatype? I don't think the former is possible. For the latter, I think the suggestions given here should suffice. – tiago Dec 16 '12 at 10:15
  • 3
    I agree with tiango, from the question I would say this is incorrect, because rounding does not guarantee to find the best integer solution (unless it is integer to begin with), but with an overdetermined system at hand, that seems unlikely. The example happens to have an exact solution though I admit. – seberg Dec 16 '12 at 11:31
  • @A.R.S. - No, as far as I know numpy doesn't have any way to solve a system of equations that's guarenteed to yield an integer solution. Internally, numpy is just calling routines from `LAPACK`, etc, so functions such as `lstsq` are inherently floating point. I think between NPE's answer to determine if a solution exists and seberg's suggestion to try a non-linear solver, you can probably put something together, though. A rounded or truncated floating point solution might give a good starting guess for the non-linear solver, if nothing else. – Joe Kington Dec 16 '12 at 15:34
1

+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
aka.nice
  • 9,100
  • 1
  • 28
  • 40
1

I needed to do this and ended up porting a PHP program written by Keith Matthews, which you can find on http://www.numbertheory.org/php/php.html, into Python. I initially used Numpy arrays but ran into problems with integer overflows so switched to Sympy matrices, which use arbitrary precision numerical representations.

The code is released on GitHub at https://github.com/tclose/Diophantine under the MIT licence, so feel free to use it and let me know if you run into any problems (sorry it isn't better documented). The master branch uses Sympy but you can access the original Numpy implementation in the 'numpy' branch, which seems to work okay for reasonably sparse systems.

If you do end up using it for a scientific publication please cite Keith's papers (and maybe add a link to the GitHub repo).

Tom Close
  • 544
  • 6
  • 12
1

My approach is to find non-ineger solution first, then scale up to integer one

from fractions import Fraction, gcd
from functools import reduce

def lcm(a, b):
    return a * b // gcd(a, b)

def common_int(*numbers):
    fractions = [Fraction(n).limit_denominator() for n in numbers[0]]
    multiple = reduce(lcm, [f.denominator for f in fractions])
    ints = [f * multiple for f in fractions]
    divisor = reduce(gcd, ints)

    return [int(n / divisor) for n in ints]

sol = np.linalg.solve(np.array([[1, 2, 3], [2, 1, 0], [2, 1, 4]]), np.array([1., 1., 1.]))  # system of equation
# [0.3333333, 0.3333333, 0.]

common_int(sol)
# [1., 1., 0.]
Louis Szeto
  • 11
  • 1
0

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.