I stopped contributing to Sage since SageMath Inc. is a for-profit Delaware company with an actual self-appointed CEO whose top priority is to "make a lot of money". The website sagemath.com steals traffic from sagemath.org to promote SageMath Inc’s products, and SageMath Inc. takes its orders from unknown investors while free contributors to Sage have no say in the matter.

I did not spend years’ worth of nights and week-end writing open-source code so that private individuals will get richer. My government does not pay me to be the free workforce of a private for-profit company either.

Linear Programming in Sage

Nathann Cohen
nathann (this round thing) cohen (the weird 'a') gmail (same round thing) com

The Linear Programs, when practically dealing with graphs are a *GREAT* tool. There are many very good solvers around, and they are now available in Sage. This means that using them only takes a couple of lines, and writing an algorithm to solve the matching problem takes about 5 minutes, which is already more than enough.

What is a Linear Program ?

A linear program is the sum of two information :

  • A linear function, called the objective, which is to be maximized or minimized ( for example 2 x + y )
  • Linear constraints on the variables ( for example, 3 x + y < 2 and 2 x + 3 y < 8 )
The solver will then try to find a solution to the system of constraints such that the objective function is optimized, and return the values of the variables.

What is a Mixed Integer Linear Program ?

It is simply a Linear Program such that some variables are forced to take integer values instead of real values. This difference becomes very important when one learns that solving a Linear Program can be done in polynomial time while solving a general Mixed Integer Linear Program is usually NP-Complete ( = it can take exponential time, according to a widely-spread belief that P is not equal to NP)

Why is Linear Programming so useful ?

Linear Programming is very useful in many Optimization and Graph-Theoretical problems because of its wide range of expression. Most of the time, a natural Linear Program can be easily written to solve a problem whose solution will be quickly computed thanks to the wealth of heuristics already contained in Linear Program Solvers. It is often hard to theoretically find out the execution time of a Linear Program, though they give very interesting results in practice.

For more information, you can consult the Wikipedia page dedicated to Linear Programming : http://en.wikipedia.org/wiki/Linear_programming

How can I solve a linear program using Sage ?

I have been trying to keep the use of MILP in Sage as easy as possible, and to have it look like what we draw on blackboard. If it takes two lines to write it on a blackboard, then it should be the same in Sage.

Sage can solve Linear Programs or Mixed Integer Linear Programs through the class 'MixedIntegerLinearProgram' defined in 'sage.numerical.mip'. To illustrate how it can be used, we will try to solve one easy problem

First, we define our MixedIntegerLinearProgram object as a maximization problem

  sage: p=MixedIntegerLinearProgram( maximization=True ) 

We then define the objective using p[1] and p[2] as variables. Not that we could have picked `p["x1"]` and `p["x2"]`, or even `p[("x",1)]` and `p[("x",2)]` : any hashable Sage object is a valid key, as it would be in an usual Python dictionary !

  sage: p.set_objective( 2*p[1]+p[2] ) 

The two constraints :

  sage: p.add_constraint( 3*p[1]+4*p[2], max = 2.5 ) 
  sage: p.add_constraint( 1.5*p[1]+0.5*p[2], min = 0.5, max = 4 ) 

Our problem is now defined, and we can now ask Sage to solve it.
( be assured that you have installed the Sage package for Coin-OR CBC or GLPK. If you do not, see the Solvers section at the bottom of this page )

  sage: p.solve()

The value returned by the ``solve`` method is the optimal value of the objective functions. In order to obtain the optimal values, you have to use the ``get_value``

  sage: print "The optimal values are x_1="+str(p.get_values(p[1]))+"x_2="str(p.get_values(p[2]))

All in all, here is how to solve this system

  sage: p=MixedIntegerLinearProgram( maximization=True )
  sage: p.set_objective( 2*p[1]+p[2] )
  sage: p.add_constraint( 3*p[1]+4*p[2], max=2.5 )
  sage: p.add_constraint( 1.5*p[1]+0.5*p[2], max=4,min=0.5 )
  sage: p.solve()
  sage: print "The optimal values are x_1 = "+str(p.get_values(p[1]))+", x_2 = "+str(p.get_values(p[2]))
  The optimal values are x_1 = 0.833333333333, x_2 = 0.0

Dictionaries of variables in ``MixedIntegerLinearProgram``

For more complex Linear Programs, it is sometimes useful to associate many parameters to the same object. Instead of using just one dictionary, it is possible to create many of them using the method ``new_variable`` : if several objects `B1, B2, ...` all have both a "value" and a "cost", the two being variables of the same linear program, you can define two dictionaries to store them

  sage: value = p.new_variable()
  sage: cost = p.new_variable()

It is now possible to define constraints mixing the two different kinds of variables

  sage: p.add_constraint( cost[ B1 ] +  3*value[ B2 ], max = 9)

Some famous examples

Vertex Cover in a graph

In the Vertex Cover problem, we are given a graph `G` and we want to find a subset `S` of its vertices of minimal cardinality such that each edge `e` is incident to at least one vertex of `S`. In order to achieve it, we define a binary variable `b_v` for each vertex `v`.

In the linear program, the syntax is exactly the same :
  sage: g=graphs.PetersenGraph() 
  sage: p=MixedIntegerLinearProgram(maximization=False) 
  sage: b=p.new_variable() 
  sage: for (u,v) in g.edges(labels=None): 
  ...:      p.add_constraint(b[u]+b[v],min=1) 
  sage: p.set_objective(sum([b[v] for v in g]))
  sage: p.set_binary(b) 
  sage: p.solve()

  sage: b = p.get_values(b)
  sage: print b
  {0: 0.0, 1: 1.0, 2: 0.0, 3: 1.0, 4: 1.0, 5: 1.0, 6: 1.0, 7: 1.0, 8: 0.0, 9: 0.0}

  sage: print [v for v in g if b[v]==1]
  [1, 3, 4, 5, 6, 7]

Maximum matching in a Graph

In the maximum matching problem, we are given a graph `G`, and we are looking for a set of edges `M` of maximum cardinality such that no two edges from `M` are adjacent :

Here is how this is solved through Sage on a Petersen Graph :
  sage: g=graphs.PetersenGraph() 
  sage: p=MixedIntegerLinearProgram() 
  sage: b=p.new_variable() 
  sage: B = lambda x, y : b[Set([x,y])]

  sage: for u in g: 
  ...:    p.add_constraint(sum([B(u,v) for v in g.neighbors(u)]),max=1) 
  sage: p.set_objective(sum([B(u,v) for (u,v) in g.edges(labels=None)]))

  sage: p.solve()
  sage: b = p.get_values(b)
  sage: print [(u,v) for (u,v) in g.edges(labels=None) if B(u,v) == 1]
  [(0, 5), (1, 6), (2, 7), (3, 8), (4, 9)]

And the next step is ``p.solve()`` !

Maximum Independent Set

A maximum independent set in a graph is a maximum set of vertices which are not connected to each other. It can be easily formulated as a LP program

... and here is how you can do it in Sage :
  sage: G = graphs.PetersenGraph()
  sage: LP = MixedIntegerLinearProgram(maximization=True)
  sage: b = LP.new_variable()

  sage: # We define the objective
  sage: LP.set_objective(sum([b[v] for v in G]))

  sage: # For any edge, we define a constraint
  sage: for (u,v) in G.edges(labels=None):
  ...       LP.add_constraint(b[u]+b[v],max=1)

  sage: # The variable b is binary
  sage: LP.set_binary(b)
  sage: # The .solve() functions returns the objective value
  sage: LP.solve()
  sage: b_sol = LP.get_values(b)
  sage: print b_sol
  {0: 1.0, 1: 0.0, 2: 1.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 1.0, 9: 1.0}


Sage solves linear programs by calling specific libraries. Two are available for the moment :

  • GLPK (website) : A Linear Program solver from GNU
  • CBC (website) : Integer Linear Program solver from COIN-OR
  • CPLEX (website) : Integer Linear Program solver from IBM ILOG. CPLEX is proprietary but free for researchers and students !

To install them if they are not available on your installation of Sage, type :

  sage: # To install GLPK 
  sage: install_package('glpk') 
  sage: # To install Coin-OR Branch and Cut ( CBC ) 
  sage: install_package('cbc') 
Depending on the version you are using, it may be a good idea to run in a console the following line
  sage -b
So that Sage notices a new package has been installed !

Get help and documentation

Short of the two very useful tips discussed in the previous section (the "Tabulation" key and the question mark "?" at the end of a line), there are three ways to obtain help :
  • The Sage-support Google group (here)
  • Sage's documentation page (here), containing for example :
  • Sage's IRC Channel : sage-support@irc.freenode.net
  • and .... do not hesitate to send me an email if you have any question ! ;-)