Ask Your Question
1

Solving linear inequality systems using Sage

asked 2022-04-29 11:23:28 +0200

klaaa gravatar image

updated 2022-04-30 17:06:19 +0200

Given a matrix A with rational entries, for example the 36x7 matrix:

A=[ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ]

Is there an easy way using Sage (or another computer algebra system) to check whether there is a solution to to the inequality systems Ax>0 (meaning every entry of the vector of Ax is strictly positive) in rational (or even real) numbers?

So the input should be the matrix A and the output should be whether there is a solution of Ax>0 or not in the rationals numbers (and if possible give a solution if there is one).

edit retag flag offensive close merge delete

Comments

sage: Polyhedron?

FrédéricC gravatar imageFrédéricC ( 2022-04-29 11:47:52 +0200 )edit

Thank you for that comment. I am not sure what you mean. Can you give an example for the much smaller list W=[ [ [ 1/2x1+1/2x2 ], [ x2 ] ], [ [ 1/2x2+1/2x3 ], [ x3 ] ], [ [ 1/3x1+1/3x2+1/3x3 ], [ 1/2x2+1/2*x3 ] ] ] using Sage?

klaaa gravatar imageklaaa ( 2022-04-29 11:51:02 +0200 )edit

1 Answer

Sort by » oldest newest most voted
1

answered 2022-04-29 16:16:18 +0200

Max Alekseyev gravatar image

updated 2022-04-30 23:12:44 +0200

To find just one solution you can pose the problem as MILP with .set_objective(None). See documentation and examples at https://doc.sagemath.org/html/en/refe...

To deal with strict inequalities, one can multiply them by LCM of the denominators (turning both sides into integers) and add 1 to the smaller side to turn the inequality into non-strict. For example,

1/9*x2+2/9*x3+2/9*x4+2/9*x5+1/9*x6+1/9*x7 > 1/4*x2+1/4*x3+1/4*x4+1/4*x5

is equivalent to

4*(x2+2*x3+2*x4+2*x5+x6+x7) >= 9*(x2+x3+x4+x5) + 1.

As for solving in rational numbers, this can be done like this:

from sage.numerical.mip import MIPSolverException

A = Matrix(QQ, [ [ 1/10, 1/10, -1/20, -1/20, -1/20, 1/10, -3/20 ], [ 0, -4/21, -1/21, -1/21, 1/7, 0, 1/7 ], [ -3/88, 5/88, 1/44, 5/88, -3/88, -3/88, -3/88 ], [ 1/13, -3/65, 2/65, -3/65, -3/65, 1/13, -3/65 ],
[ -5/36, -1/36, -1/36, -1/36, 1/9, 0, 1/9 ], [ 1/6, 1/6, -1/6, 1/6, 0, 0, -1/3 ], [ 0, -2/3, 1/3, 0, 0, 0, 1/3 ], [ -1/2, 1/2, 0, 0, 0, 0, 0 ],
[ 0, 0, -1/12, -1/12, -1/12, 1/4, 0 ], [ 0, -5/36, -1/36, -1/36, -1/36, 1/9, 1/9 ], [ 1/13, 1/91, -5/91, -5/91, 1/91, 1/13, -6/91 ], [ 1/15, -1/30, -1/15, 1/30, -1/30, 1/15, -1/30 ],
[ -3/208, 7/208, 1/52, 7/208, -3/104, -3/208, -3/104 ], [ 1/153, -7/153, 2/153, -7/153, 1/153, 1/17, 1/153 ], [ 1/11, -3/44, 1/44, -3/44, 1/11, 0, -3/44 ], [ 1/7, -1/21, -1/21, 1/7, 0, 0, -4/21 ],
[ 1/4, 1/4, 1/4, 0, 0, 0, -3/4 ], [ 0, 0, 0, -1/6, -1/6, 1/3, 0 ], [ 0, 0, -1/6, -1/6, 1/3, 0, 0 ], [ 0, -1/56, -1/28, 3/28, -1/56, -1/56, -1/56 ],
 [ -1/132, -1/66, 3/44, -1/66, -1/66, -1/132, -1/132 ], [ -3/56, 1/56, -1/28, -1/28, 1/56, 1/14, 1/56 ], [ 1/10, 0, -1/10, 0, 1/10, 0, -1/10 ], [ -1/132, -1/66, -1/44, -1/66, -1/132, 1/12, -1/66 ],
 [ -1/24, 1/12, 1/12, -1/24, -1/24, 0, -1/24 ], [ -1/20, -1/20, -1/20, 1/5, 0, 0, -1/20 ], [ 0, 0, 0, 0, -1/2, 1/2, 0 ], [ 0, 0, 0, -1/2, 1/2, 0, 0 ],
 [ 0, 0, -1/2, 1/2, 0, 0, 0 ], [ 0, -1/42, 5/42, -1/42, -1/42, -1/42, -1/42 ], [ -6/55, -1/55, -1/55, -1/55, -1/55, 1/11, 1/11 ], [ 1/8, 1/8, -1/12, -1/12, 1/8, 0, -5/24 ],
 [ 0, -3/10, -1/10, 1/5, 0, 0, 1/5 ], [ -2/63, 5/63, 5/63, -2/63, -2/63, -2/63, -2/63 ], [ -1/30, -1/30, -1/30, -1/30, 1/6, 0, -1/30 ], [ -1/42, -1/42, -1/42, -1/42, -1/42, 1/7, -1/42 ] ])
#print(A)

eps = 1e-6

mylp = MixedIntegerLinearProgram( solver='ppl' )
mylp.set_objective( None )
x = mylp.new_variable()
for r in A.rows():
    mylp.add_constraint( sum(r[i]*x[i] for i in range(len(r))) >= eps )

try:
    mylp.solve()
    print('Solution:', mylp.get_values(x))
except MIPSolverException:
    print('No solution')

Here we introduced a small positive gap eps to turn strict inequalities into non-strict. However, we get no solution even if further decrease eps. This suggests that the original problem has no solution. To verify it we can construct the LP problem for eps = 0 as above and check the volume of the corresponding polyhedron:

mylp.polyhedron(base_ring=QQ).volume()

Since the volume is zero, the polyhedron defined by $Ax\geq 0$ has no internal points and thus the inquality $Ax>0$ has no solutions.

edit flag offensive delete link more

Comments

Thank you very much. Is there by the way a theoretical result that decides whether a system of inequalities of the form Ax>0 for a matrix with integer entries A has a solution?

klaaa gravatar imageklaaa ( 2022-04-29 17:58:47 +0200 )edit
1

Integer programming in general is NP-hard - see https://en.wikipedia.org/wiki/Integer...

Max Alekseyev gravatar imageMax Alekseyev ( 2022-04-29 18:32:51 +0200 )edit

Thanks, that is interesting. But is there some way to find out what the best possible algorithm known at the moment is to check whether the problem Ax>0 has a solution (in rational numbers or integers or even real numbers) for a given matrix A with rational entries? I looked on some books on optimation but they are a bit old.

klaaa gravatar imageklaaa ( 2022-04-30 12:41:09 +0200 )edit

Solving in rational or real numbers is much simpler than solving in integers. In this case we talk about linear programming(LP) rather than ILP. Sage provides access to a variety of (I)LP solvers, both freeware and commercial, which can can solve your LP or may be even ILP problem efficiently. Check out this tutorial: https://doc.sagemath.org/html/en/them...

Max Alekseyev gravatar imageMax Alekseyev ( 2022-04-30 15:44:52 +0200 )edit

Thank you very much again! The question in the link still looks a bit different when searching for rational solutions. I simplified my question a bit. Can you give a concrete example how to do it with the matrix A as in the question?

klaaa gravatar imageklaaa ( 2022-04-30 17:07:31 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2022-04-29 11:23:28 +0200

Seen: 150 times

Last updated: Apr 30