pplp / blog / blog.html

Full commit
<H1>Playing with Linear Programming on PyPy</H1>
Fancy hi-level interfaces often come with a high runtime overhead
making them slow. Here is an experiment with building such an
interface using constructions that PyPy should be good at
optimizing. The idea is to allow the JIT in PyPy to remove the
overhead introduced by using a fancy high-level python interface
on top of a low-level C interface. The application considered is
<a href="">Linear
programming</a>. It is a tool used to solve linear optimization
problems. It can for example be used to find the nonnegative values
x, y and z that gives the maximum value of

<p />
<center><div class="eq"><img alt="" src="eqsource1.png" /></div></center>

<p />
without violating the constraints

<p />
<center><div class="eq"><img alt="" src="eqsource2.png" /></div></center>
<center><div class="eq"><img alt="" src="eqsource3.png" /></div></center>
<center><div class="eq"><img alt="" src="eqsource4.png" /></div></center>

<p />
There exists general purpose solvers for these kind of problems that
are very fast and can literally handle millions of variables. To use
them however the problem has to be transformed into some specific
matrix form, and the coefficients of all the matrices
has to be passed to the solver using some API. This transformation is
a tedious and error prone step that forces you to work with matrix
indexes instead of readable variable names. Also it makes maintaining
an implementation hard since any modification has to be transformed

<p />
The example above comes from the manual of
the <a href="">glpk</a> library. That
manual continues by describing how to convert this problem into the
standard form of glpk (which involves introducing three new variables)
and then gives the <a href="sample.c">c-code</a> needed to call the
library. Relating that c-code to the problem above without the
intermediate explanation of the manual is not easy. A common
solution here is to build a hi-level interface that allows a more
natural way of defining the matrices and/or allow the equations to be
entered symbolically. Unfortunately, such interfaces often become
slow. For the benchmark below for example, 
<a href="">cvxopt</a>
requires 20 minutes to setup a problem that takes 9.43 seconds to solve
(this seems a bit extreme, am I doing something wrong?).

<p />
The high-level interface I constructed on top of the
<a href="">glpk</a> library is 
<a href="">pplp</a> and it allows
the equations to be entered symbolically. The above problem can be
solved using
    lp = LinearProgram()
    x, y, z = lp.IntVar(), lp.IntVar(), lp.IntVar()
    lp.objective = 10*x + 6*y + 4*z
    lp.add_constraint( x + y + z &lt;= 100 )
    lp.add_constraint( 10*x + 4*y + 5*z &lt;= 600 )
    lp.add_constraint( 2*x + 2*y + 6*z &lt;= 300 )
    lp.add_constraint( x &gt;= 0 )
    lp.add_constraint( y &gt;= 0 )
    lp.add_constraint( z &gt;= 0 )

    maxval = lp.maximize()
    print maxval
    print x.value, y.value, z.value

<p />
To benchmark the API I used it to solve a 
<a href="">minimum-cost
  flow problem</a> with 154072 nodes and 390334 arcs. The C library
  needs 9.43 s to solve this and the pplp interface adds another 5.89
  s under PyPy and 28.17 s under CPython. A large amount of time is
  still spend setting up the problem, but it's a significant
  improvement over the 20 minutes required on CPython by
  <a href="">cvxopt</a>. It is
  probably not designed to be fast on this kind of benchmark. I have
  not been able to get cvxopt to work under PyPy. The benchmark used is
  available <a href="">here</a>