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
On Friday, September 14, 2012 5:18:15 PM UTC+2, 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.
>
> Victor
>
> On Friday, September 14, 2012 5:46:01 AM UTC-4, Nathann Cohen wrote:
>>
>> 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
>>
>
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To post to this group, send email to sage-support@googlegroups.com.
To unsubscribe from this group, send email to
sage-support+unsubscr...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.
