Goal

This project provides several algorithms implemented in python to solve linear programs of the form

where A_e and A_i are sparse matrices

The different algorithms that are implemented here are documented in the pdf:

• a dual coordinate ascent method with exact line search
• a dual gradient ascent with exact line search
• a first order primal-dual algorithm adapted from chambolle pock [2]
• three methods based on the Alternating Direction Method of Multipliers [3]

Note These methods are not meant to be efficient methods to solve generic linear programs. They are simple and quite naive methods I implemented while exploring different possibilities to solve very large sparse linear programs that are too big to be solved using the standard simplex method or standard interior point methods.

This project also provides:

• a python implementation of Mehrotra's Predictor-Corrector Pimal-Dual Interior Point method.
• a python class SparseLP (in SparseLP.py) that makes it easier to build linear programs from python
• methods to convert between the different common forms of linear programs (slack form, standard form etc),
• methods to import and export the linear program from and to standard file formats (MPS), It is used here to run netlib LP problems. Using mps files, one can upload and solve LP on the neos servers.
• a simple constraint propagation method with back-tracking to find feasible integer values solutions (for integer programs)
• an interface to the OSQP solver [11]
• interfaces to other solvers (SCS, ECOS, CVXOPT) through CVXPY
• interfaces to other LP and MILP solvers (CLP, CBC, MIPLC, GLPSOL, QSOPT) using mps text files

Build and test status

Installation

using pip

sudo pip install git+git://github.com/martinResearch/PySparseLP.git

otherwise you can dowload it, decompress it and compile it locally using

python setup.py build_ext --inplace

If you want to be able to run external solvers using mps files in windows then download the following executables and copy them in the solvers\windows subfolder

• clp.exe
• cbc.exe
• mps_mipcl.exe
• glpsol.exe and glpk_4_63.dll
• qssolver.exe

LP problem modeling

This library provides a python class SparseLP (in SparseLP.py) that aims at making it easier to build linear programs from python. It is easy to derive a specialize class from it and add specialized constraints creations methods (see potts penalization in example 1). SparseLP is written in python and relies on scipy sparse matrices and numpy matrices to represent constraint internally and for its interface. There is no variables class binding to c++ objects. This makes it potentially easier to interface with the python scientific stack.

Debuging

Constructing a LP problem is often error prone. If we can generate a valid solution before constructing the LP we can check that the constraints are not violated while we add them to the LP using the method check_solution. This make it easier to pin down which constraint is causing problem. We could add a debug flag so that this check is automatic done as we add constraints.

Other modeling tools

Other libraries provide linear program modeling tools:

• CyLP. It uses operator overloading so that we can use notations that are close to the mathematical notations. Variables are defined as 1D vectors wich is sometimes not convenient.
• GLOP. The variables are added as scalars, one at a time, instead of using arrays, which make the creation of large LPs very slow in python.
• PuLP. The variables are added as scalars, one at a time, instead of using arrays, which make the creation of large LPs very slow in python.
• Pyomo
• CVXPY

The approach taken here is lower level than these tools (no variable class and no operator overloading to define the constraints) but provides more control and flexibility for the user to define the constraints and the objective function. It is made easy by using numpy arrays to store variables indices.

Examples

Image segmentation

We consider the image segmentation problem with Potts regularization:

with E the list of indices of pairs of neighbouring pixels and c a cost vector that is obtained from color distribution models of the two regions. This problem can be rewritten as a linear program by adding an auxiliary variable d_ij for each edge with the constraints

 

This problem can be more efficiently solved using graph-cuts than with a generic linear program solver but it is still interesting to compare performance of the different generic LP solvers on this problem.

from pysparselp.example_pott_segmentation import run
run()

Here are the resulting segmentations obtained with the various LP solvers, with the same random data term with the optimizations limited to 15 seconds for each solver. convergence curves

Note that instead of using a simple Potts model we could try to solve the LP from [5]

Sparse inverse convariance matrix

The Sparse Inverse Covariance Estimation problem aims to find a sparse matrix B that approximate the inverse of Covariance matrix A.

Let denote f the fonction that take a matrix as an input an yield the vector of coefficient of the matrix in row-major order. Let b=f(B) we have f(AB)=Mb with M=kron(A, I_d) The problem rewrites

We take inspiration from this scikit-learn example here to generate samples of a gaussian with a sparse inverse covariance (precision) matrix. From the sample we compute the empirical covariance A and the we estimate a sparse inverse covariance (precision) matrix B from that empirical covariance using the LP formulation above.

from pysparselp.example_sparse_inv_covariance import run
run()

L1 regularised multi-class SVM

Given n examples of vector-class pairs (x_i,y_i), with x_i a vector of size m and y_i an integer representing the class, we aim at estimating a matrix W of size k by m that allows to discriminate the right class, with k the number of classes. We assume that the last component of x_i is a one in order to represent the offset constants in W. we denote W_k the k^th line of the matrix W

By adding auxiliary variables in a matrix S of the same size as the matrix W we can rewrite the absolute value as follow:

We obtain the LP formulation:

The example can be executed using the following line in python

from pysparselp.example_sparse_l1_svm import run
run()

The support vectors are represented by black circles.

Bipartite matching

Bipartite matching can be reformulated as an integer linear program:

We relax it into an continuous variables LP.

from pysparselp.example_bipartite_matching import run
run()

K-medoids

Given n points we want to cluster them into k set by minimizing

with d_ij the distance between point i and point j This can be reformulated as an integer program:

 

We relax it into a continuous variables LP using  

from pysparselp.example_kmedians import run
run()