2nd try: Integer solutions for systems of linear equations

[moritz]

Note: The first time I tried to post this it didn't show on the google groups UI, so I try again..

-----------

Let's consider a system of linear equations with integer coefficients:

n=3

m=5

A=random_matrix(ZZ,n,m)

b=random_vector(ZZ,n)

What is the best way to find an integer vector x such that Ax=b?

A rational solution is of course found by A.solve_right(b). 

A random example:

 

A=matrix(ZZ,3,5,[-1,  0, -1, -1,  1, 0, -2, -1,  1, -8,-1,  0,  2, -8, -2])

b=vector(ZZ,3,[-1,-1,15])

A.solve_right(b) gives (-13/3, -13/6, 16/3, 0, 0), although there are integer solutions,

e.g.: (-221, 13, 144, 65, -13)

Is there also a good way to get all integer solutions? In the example above that would be something like:

[-10*t0 - 30*t1 - 221, -2*t0 + 3*t1 + 13, 7*t0 + 19*t1 + 144, 3*t0 + 9*t1 + 65, -2*t1 - 13]

for integer variables t0 and t1.

The best way I could think of is to implement the algorithm from section 4.5.2 from Knuth's TAOCP. (This is also where the solutions above come from)

Is there a simpler way to do that in sage?

Of course to we could formulate the same question in terms of linear equations instead of matrices, but then solve() also gives rational, non-integer solutions. 

Maybe the right thing to use would be isolve from maple, but the interface is broke ( see http://trac.sagemath.org/sage_trac/ticket/12295 ).

To consider the problem as linear program and use MixedIntegerLinearProgram()  with integer constrains works, but it is very very slow for larger systems.

Any help will be appreciated,

    moritz

[P Purkayastha]

On 09/13/2012 03:46 PM, moritz wrote:
> Note: The first time I tried to post this it didn't show on the google
> groups UI, so I try again..

I can't help you with your math problem, but I did receive your earlier
message via gmane, although it doesn't show up in google groups for me too.

http://thread.gmane.org/gmane.comp.mathematics.sage.support/29489

[luisfe]

On Thursday, September 13, 2012 9:46:38 AM UTC+2, moritz wrote:

Note: The first time I tried to post this it didn't show on the google groups UI, so I try again..

-----------

Let's consider a system of linear equations with integer coefficients:

What is the best way to find an integer vector x such that Ax=b?

A rational solution is of course found by A.solve_right(b). 

Look for information about the Hermite normal form of an integer matrix.

[Nathann Cohen]

Helloooooooooooooo !!!

To consider the problem as linear program and use MixedIntegerLinearProgram()  with integer constrains works, but it is very very slow for larger systems.

Well, your problem is *precisely* what Integer Linear Program solvers are written for, so I guess that using them is the good way to go unless you plan on using some properties of the matrices you generate (and that the LP solvers would not notice) to solve your equation.

They are indeed slow in some cases, but we have to work with the tools available, or rather those we know about. ILP's what I would try to do myself -- if you find some other way to solve these equations please let me know, for I use them very often and I would be delighted to solve my graph problems faster too :-)

By the way, in Sage MixedIntegerLinearProgram uses GLPK by default to solve these problems, but you can also ask it to solve them with Gurobi or CPLEX, which are usually *MUCH* faster.

Good luuuuuuuuuck ! :-)

Nathann

[Victor Miller]

Since the OP's problem has no inequalities (such as requiring that all integers in question are non-negative), it is solved by using Hermite normal form.

If A is an m by n integer matrix, the Hermite normal form of A is an upper triangular integer matrix H (also m by n), along with an m by m integer matrix of determinant 1 (i.e. it's invertible) so that H = U A.  So the original equation A x = b becomes (after multiplying by U): H x = U b, for which it's easy to get the general solution (just back solve.  If you ever get a non-integer by dividing there are no solutions).  Since U is invertible, multiplying the equation by U doesn't introduce spurious solutions.

Victor

[Nathann Cohen]

> Since the OP's problem has no inequalities (such as requiring that all
> integers in question are non-negative), it is solved by using Hermite normal
> form.
>
> If A is an m by n integer matrix, the Hermite normal form of A is an upper
> triangular integer matrix H (also m by n), along with an m by m integer
> matrix of determinant 1 (i.e. it's invertible) so that H = U A.  So the
> original equation A x = b becomes (after multiplying by U): H x = U b, for
> which it's easy to get the general solution (just back solve.  If you ever
> get a non-integer by dividing there are no solutions).  Since U is
> invertible, multiplying the equation by U doesn't introduce spurious
> solutions.

**Applause**

Well... Thank you for teaching me that :-)

Nathann

[Jason Grout]

On 9/14/12 10:18 AM, Victor Miller wrote:
> Since the OP's problem has no inequalities (such as requiring that all
> integers in question are non-negative), it is solved by using Hermite
> normal form.
>
> If A is an m by n integer matrix, the Hermite normal form of A is an
> upper triangular integer matrix H (also m by n), along with an m by m
> integer matrix of determinant 1 (i.e. it's invertible) so that H = U A.
> So the original equation A x = b becomes (after multiplying by U): H x =
> U b, for which it's easy to get the general solution (just back solve.
> If you ever get a non-integer by dividing there are no solutions).
> Since U is invertible, multiplying the equation by U doesn't introduce
> spurious solutions.

And for completeness, you can get hermite form by using the
.hermite_form method:


sage: a=random_matrix(ZZ,5)
sage: a.hermite_form()
[ 1 0 0 0 35]
[ 0 1 0 0 73]
[ 0 0 1 0 40]
[ 0 0 0 1 51]
[ 0 0 0 0 103]
sage: a.hermite_form(transformation=True)
(
[ 1 0 0 0 35] [ -37 -108 28 -10 97]
[ 0 1 0 0 73] [ -77 -224 58 -21 201]
[ 0 0 1 0 40] [ -44 -128 33 -12 115]
[ 0 0 0 1 51] [ -59 -172 44 -16 154]
[ 0 0 0 0 103], [-110 -321 83 -30 288]
)


Thanks,

Jason

[moritz]

Thanks for the reference to Hermite normal forms. I haven't quite figured out how to use it to solve the problem. Either I misunderstood something or it is not quite as easy as you describe it. Take this example

A=matrix(([1,5,0],[0,6,4]))

b=vector([3,1])

c=vector([5,2])

and consider Ax=b and Ax=c

so A is already in Hnf. But simple back solving does not do the trick: For Ax=b in the last row the equation becomes

0*x[0]+6*x[1]+4*x[2]==1

and according to your method, we get no solution (which is true), since we get a non-integer by dividing by 6.

For Ax=c, the last row reads:

0*x[0]+6*x[1]+4*x[2]==2

and again we could think that there is no solution since 2/6 is not an integer, but if we allow x[2] to be non-zero, we get the solution x[1] = 1, x[2] = -1.

(This can't be done for Ax=b, since the equation has no solution mod 2)

In generally, when "back solving" the equations maybe one has to calculate the gcd of the rows, and then find an integer solution, using some kind of extended Euclid's algorithm, if the gcd is 1. (Then a solution exists according to the Chinese remainder theorem..)

Anyhow I found a different solution: the smith normal form! (Which is already implemented in sage..)

def smithsolve(A,b):

    D,U,V=A.smith_form()

    c=U*b

    d=D.diagonal()

    y=vector(ZZ,D.ncols())

    for i in [len(d),..,D.nrows()-1]:

        if c[i]:

            return 'no integer solution'

    for i in range(len(d)):

        if d[i]==0:

            if c[i]==0:

                y[i]=0

            else:

                return 'no integer solution'

        else:

            q=c[i]/d[i]

            if q in ZZ:

                y[i]=q

            else:

                return 'no integer solution'

    x=V*y

    return x

n=3

m=6

A=random_matrix(ZZ,n,m)

b=random_vector(ZZ,n)

x=smithsolve(A,b)

print A,b, x, A*x==b

Maybe an algorithm that uses the Hermite normal form would be faster.

One questions remains:

Should sage provide the functionality to give integer solutions of linear Diophantine equations if one asked for it, e.g. in the solve_right command? Wouldn't it be nice if the solve_right command could take an extra argument like "base ring" and try to solve it in that ring. For example something along the lines:

A.solve_right(b,base_ring=ZZ)

moritz