[SciPy-Dev] [SciPy-User] Linear Programming via Simplex Algorithm

[Christopher Jordan-Squire]

Rob--Thanks for posting the code. I've been hoping an LP solver would
make it into scipy.optimize for awhile. Like Rob said, they're
important both on their own and as subproblems in other algorithms,
such as integer LP arising from approximation algorithms for
combinatorial optimization problems. Successive linear programming has
also had a lot of success in some industries. Rob's also right that
it's generally better to use an LP solver on an LP problem than a
non-linear solver. If comparing implementations of similar quality,
then an LP solver will usually be faster and more accurate.

Anyways, I wouldn't consider myself an expert, but I am semi-informed.
For people not familiar with linear programming, it's probably good to
mention a few important facts about LP's.

The tl;dr is LP solvers are really complicated. Including a basic one
in Scipy is possible, especially with iterative improvements over
time, but having even a basic LP solver with all the major options
included for a standard LP solver would be a pretty huge undertaking.
Details about why are below. Having just the simplex (and dual
simplex? What happened to that PR?) would be nice, but it's important
to be up front about them being relatively untested, almost certainly
not very robust, and (I assume) using dense matrices instead of sparse
matrices.

<general LP info>

--Existing open source implementations are generally about an order of
magnitude slower than commercial implementations. [1] [2] show this,
but those are specifically talking about integer and mixed integer
LP's.

--The main commercial implementations are Cplex, Gurobi, Mosek, and
Xpress. As far as I know their interfaces are rather dissimilar,
and--unfortunately--many OR people use a bafflingly large array of
DSL's to build their LPs, such as AMPL, GAMS, Mosle (Xpress), and OPL
(Cplex).

--There are 3 major algorithms for solving general LP's: simplex, dual
simplex, and barrier. The simplex and dual simplex methods will, for
large-scale problems, use sparse direct linear algebra solvers, while
barrier methods can, depending on the implementation, use sparse
direct or iterative solvers.

--There are also a whole host of algorithms for efficiently solving
LP's with more structure, such as network flow LP's. These include
primal-dual method (different from primal-dual interior point method)
and auction algorithms. As I understand it, commercial LP solvers
generally do some checking to recognize if that kind of structure is
present.

--Mature LP implementations include support for large, sparse LP's,
which uses some LP-specific sparse matrix techniques [3]. They also
include presolving, but the details of that are usually proprietary.
These are a huge deal, but a standard reference says "Relatively
little information about presolving techniques has appeared in the
literature, in part because they have commercial value as an important
component of linear programming software." [4] Similarly, any time
iterative solvers are used the choice of preconditioner can be a big
deal.

--Mature implementations are also good at dealing with various kinds
of degeneracy and detecting infeasibility and unboundedness.

--All of the above refers only to LP's. Integer LP's and Mixed Integer
LP's get (much) more complicated. For example, Frank Hutter at UBC has
several papers on using an optimization algorithm just to find the
best set of configurations for Cplex for a mixed integer linear
program.

[1]http://www.gurobi.com/pdf/Benchmarks.pdf, slide 16
[2]http://scip.zib.de/
[3]http://www.stanford.edu/class/msande318/notes/notes05-updates.pdf
[4]Nocedal and Wright, Nonlinear Optimization 2nd edition, page 388

</general LP info>

Here are some specific comments on Rob's code.

--What is the cycle variable?
--I'm not a fan of matlab, but if you're going to call it linprog it
might be nice to have the same interface at matlab's
--I can't think of a good reason, either theoretical or practical, to
use both A_lb and A_up
--Similarly, I think the objective type (min or max) isn't worth
including since you can just take the negative of the cost vector
--linprog looks like it contains a lot of code that looks very
similar. Could that be refactored out into some loops and/or function
calls, to increase clarity, conciseness, and maintainability?
--What is an artificial variable?
--What is it you're doing to prevent cycling? (i.e., can you provide a
reference? I believe there are a number of different cycling rules, so
there's some ambiguity.)

 -Chris

On Wed, Oct 23, 2013 at 10:38 AM, Pauli Virtanen <pav at iki.fi> wrote:
> Hi,
>
> 23.10.2013 04:35, Rob Falck kirjoitti:
> [clip]
>> I've spent some time recently polishing a simplex-based linear
>> programming function.  I've seen traffic from time to time about including
>> such functionality in scipy.optimize but it always seems to have been
>> closed without inclusion.
>
> One important question: is this algorithm regarded as useful and robust
> enough by people in the know?
>
> Are there existing more sophisticated LP solver codes with compatible
> licences?
>
> Does the LP simplex method win over nonlinear solvers such as COBYLA?
>
> scipy.optimize currently contains many "naive" solvers, which is not a
> completely happy situation. Of course, something is better than nothing,
> but if possible, I'd prefer one sophisticated code over many naive
> methods. I'm not an expert in LP, so I can't judge very well myself here.
>
> --
> Pauli Virtanen
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev