[SciPy-Dev] Interested in contributing BILP solver and interior-point method for linprog; SciPy-dev search down?

[Matt Haberland]

Thanks for pointing out the more restrictive license of cvxopt, Gaël. Along
those lines:
lp_solve <http://lpsolve.sourceforge.net/5.5/LGPL.htm> is Lesser GPL.
GLPK <https://en.wikipedia.org/wiki/GNU_Linear_Programming_Kit> is GPLv3
(as one would expect)
COIN-OR <https://www.coin-or.org/Clp/faq.html> is CPL.
QSOPT <http://www.math.uwaterloo.ca/~bico/qsopt/ex/> is GPLv2.1+ or for
research use one (user's choice).
The rest of the codes are commercial, I believe.

As mentioned in my last email, I wanted to compare the performance of
linprog and the interior point code, which I'll refer to as linprog_ip. The
benchmark problems used were scalable and typically randomly-generated.
They were selected purely for their availability/ease of implementation
(except for the Klee-Minty Hypercube, which is the bane of simplex solvers.
It is possible that many of these tests are particularly challenging to
simplex solvers, but the only one I know to punish simplex is the
Klee-Minty hypercube.) I'd love to test with the more practical NETLIB
<http://www.cuter.rl.ac.uk/Problems/netlib.shtml> benchmarks but these are
in MPS format and I haven't looked into how to convert yet.

The results are available at:
https://www.dropbox.com/s/8vfye61kvn9h3zg/linprog_benchmarks.pdf?dl=0
and details about the tests are in the postscript.

In summary, linprog_ip is significantly faster for all but the smallest
problems (those that are solved in hundredths of a second), and is more
reliable at finding solutions when they exist.. As problem size increases,
the performance gap tends to increase. For sparse problems, the ability of
linprog_ip to leverage scipy.sparse linear algebra yields great
improvement. For certain types of problems, linprog has difficulty finding
a starting point, and takes longer to fail than linprog_ip does to return
the optimal solution.

Besides the test results, please also see all open linprog bug reports: 5400
<https://github.com/scipy/scipy/issues/5400>, 6139
<https://github.com/scipy/scipy/issues/6139>, 6690
<https://github.com/scipy/scipy/issues/6690>, 7044
<https://github.com/scipy/scipy/issues/7044>. linprog_ip suffers from none
of these. Of course, it also passes all the linprog unit tests and many
more I have written.

I attribute the success of linprog_ip to the algorithm it implements and
the suitability of the algorithm to NumPy/SciPy. For instance, for lpgen_2D
with 60 constraints and 60 variables, linprog reports about 140 iterations
whereas linprog_ip finishes in 9. Each of linprog's iterations has loops
that do one elementary row operations per loop iteration, whereas at the
heart of a linprog_ip iteration is the Cholesky factorization of a matrix
and back/forward substitution for four different right hand sides (all
LAPACK via scipy.linalg.solve). Consequently, linprog_ip spends larger
chunks of time in compiled code with far fewer loops/function calls in
Python.

Please let me know if you'd like to see the interior point code included in
SciPy, and I'll integrate with linprog and figure out the submission
process. (If you are familiar with compiling from source on Windows, please
shoot me an email. After pushing through a few roadblocks, I'm stuck with a
very vague error message.)

Matt


-----


All tests were performed on a Dell XPS 13, IntelCore i5-4200 CPU at 1.6-2.3
GHz, 8GB RAM, and Windows 10 64-bit. Tests were timed using line_profiler (
kernprof). For randomized problems, each of N problems of arbitrary size (m
constraints, n variables) were generated and solved by linprog and
linprog_ip in either (random) order. For non-random problems, the same
problem of a given dimension was solved by linprog and linprog_ip in either
(random) order N times. Unless otherwise noted, the consistency of the
solutions returned by linprog and linprog_ip was checked by comparing the
objective function value using numpy.allclose with default tolerances.
Times reported are the average per hit (call to linprog or linprog_ip ) in
units of 4.46E-07 s and also as normalized by the faster of the two times
(1 is better; the other number indicates how many times longer it took).

1. lpgen_2D <https://gist.github.com/denis-bz/8647461>, the randomized test
case currently part of the linprog unit tests.
*Note that in these tests there are actually m*n variables and m+n
constraints. For example, in the last test there are 800 constraints and
160000 variables.*
(Unfortunately, I realized after testing that I left np.random.seed(0) in
the code, and thus the same problem was solved N times per test. However,
unreported testing after fixing this yielded similar results.)

2. zionts1
<http://onlinelibrary.wiley.com/doi/10.1111/1475-3995.00392/abstract>, a
randomized test with inequality constraints originally used by Zionts to
compare the performance of his criss-cross algorithm to the simplex.
Occasionally, linprog fails to find a feasible initial point, even where a
solution exist, as reported in Issue 6139
<https://github.com/scipy/scipy/issues/6139>.

3. zionts2
<http://onlinelibrary.wiley.com/doi/10.1111/1475-3995.00392/abstract>, a
randomized test with both equality and inequality constraints originally
used by Zionts to compare the performance of his criss-cross algorithm to
the simplex. For problems of moderate size, linprog fails to find a
feasible initial point, even where a solution exist, as reported in Issue
6139 <https://github.com/scipy/scipy/issues/6139>.

4. bonates
<http://onlinelibrary.wiley.com/doi/10.1111/1475-3995.00392/abstract>, a
sparse, randomized test used in "Performance evaluation of a family of
criss–cross algorithms for linear programming". Here, percentage of
nonzeros (density) is reported as p (the sparsity is 1-p). The solutions
are typically infeasible or unbounded; the time reported is that required
for the solver to reach this conclusion. linprog always terminates when it
fails to find a feasible initial point. linprog_ip reports whether the
problem is infeasible or unbounded.

5. magic_square, an LP-relaxed version of a binary integer linear program
designed to find a magic square <https://en.wikipedia.org/wiki/Magic_square>
of size D. The cost vector is random; any feasible solution represents a
magic square. The problem includes equality constraints and simple bounds
(0, 1). As formulated, the problem's equality constraint matrix is rank
deficient (in a non-obvious way), causing linprog to report success despite
the fact that the returned solution is infeasible and the objective value
is inconsistent with the solution, as reported in Issue 6690
<https://github.com/scipy/scipy/issues/6690>. To enable performance
comparison, I've lent linprog the use of linprog_ip's presolve routine,
which removes redundant rows such that linprog is able to solve the problem
correctly.

6. klee_minty <https://en.wikipedia.org/wiki/Klee%E2%80%93Minty_cube>, a
problem that elicits the worst-case performance of Dantzig's original
simplex method. Here, Bland's rule is used by linprog. (Its runtime is far
worse with the default pivoting rule.)

On Tue, Feb 21, 2017 at 11:07 PM, Gael Varoquaux <
gael.varoquaux at normalesup.org> wrote:

> > On Tue, Feb 14, 2017 at 1:15 AM, Ralf Gommers <ralf.gommers at gmail.com>
> wrote:
> >     There's an argument to be made for having a consistent set of
> well-known
> >     algorithms in SciPy, but knowing the boundary of the scope here is
> helpful
> >     especially if there are more powerful solutions like cvxopt that
> could
> >     become easy to install soon.
>
> cvxopt is GPLv3 and scipy is BSD, so cvxopt cannot be a supplement of
> scipy for everybody.
>
> Gaël
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-dev
>



-- 
Matt Haberland
Assistant Adjunct Professor in the Program in Computing
Department of Mathematics
7620E Math Sciences Building, UCLA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20170222/b3430ac7/attachment.html>