SciPy 1.0 linear programming

The simplex solver in SciPy has been somewhat of a problem. It is a text-book LP solver (tablaux based), so only suited for very small problems. Even on those small problems it is just not very reliable.

Here is a small example [1]:

Min -x2 -x3
x0 + 0.33*x2 + 0.67*x3 = 0.5
x1 + 0.67*x2 + 0.33*x3 = 0.5
x0 + x1 + x2 + x3 = 1.0

We can encode this problem as follows:

from scipy.optimize import linprog
a_eq = [[1.0, 0.0, 0.33, 0.67],
        [0.0, 1.0, 0.67, 0.33], 
        [1, 1, 1, 1]]

b_eq = [0.5, 0.5, 1.0]

c = [0, 0, -1.0, -1.0]

res1 = linprog(c=c, A_eq=a_eq, b_eq=b_eq, method='simplex')
print(res1)

The simplex solver has troubles with this.  It returns a sub-optimal solution with objective function value 0:

     fun: -0.0
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([], dtype=float64)
  status: 0
 success: True
       x: array([ 0.5,  0.5,  0. ,  0. ])
     con: array([ -9.83213511e-13,  -9.83102488e-13,  -1.96620498e-12])

SciPy 1,0 has a new interior point solver [3] which seems a bit more robust:

res2 = linprog(c=c, A_eq=a_eq, b_eq=b_eq, method='interior-point')
print(res2)

     fun: -1.000000000000884
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([], dtype=float64)
  status: 0
 success: True
       x: array([  5.41163802e-13,   5.41163802e-13,   5.00000000e-01,
         5.00000000e-01])

We see an optimal objective of -1. It also gives a warning:

D:\Anaconda3\lib\site-packages\scipy\optimize\_linprog_ip.py:623: OptimizeWarning: A_eq does not appear to be of full row rank. To improve performance, check the problem formulation for redundant equality constraints.
  warn(redundancy_warning, OptimizeWarning)

There is no good reason a Simplex solver should give incorrect results in case the A matrix is not full-rank. Usually solvers add a slack (often not physically) to each constraint, so A will never be rank deficient. The basis matrix can always be chosen such that it is both optimal and non-singular. Here is how an optimal basis (with objective -1) looks like for this problem provided by a different simplex solver:

                    LOWER          LEVEL          UPPER         MARGINAL

---- EQU e1          0.5000         0.5000         0.5000         EPS       'Non-basic'     
---- EQU e2          0.5000         0.5000         0.5000          .        'Basic'  
---- EQU e3          1.0000         1.0000         1.0000        -1.0000    'Non-basic'  

                    LOWER          LEVEL          UPPER         MARGINAL

---- VAR x0           .              .            +INF            1.0000    'Non-basic'     
---- VAR x1           .              .            +INF            1.0000    'Non-basic'  
---- VAR x2           .             0.5000        +INF             .        'Basic'  
---- VAR x3           .             0.5000        +INF             .        'Basic'  

We see equation 2 (or its slack) is made basic although the value of the slack is zero. This solution is actually both primal and dual degenerate. [4] A different solver gives:

                   LOWER          LEVEL          UPPER         MARGINAL

---- EQU e1         0.5000         0.5000         0.5000        -1.0000    'Non-basic'   
---- EQU e2         0.5000         0.5000         0.5000        -1.0000    'Non-basic'   
---- EQU e3         1.0000         1.0000         1.0000          .        'Basic'   

                   LOWER          LEVEL          UPPER         MARGINAL

---- VAR x0          .              .            +INF            1.0000    'Non-basic'    
---- VAR x1          .              .            +INF            1.0000    'Non-basic'   
---- VAR x2          .             0.5000        +INF             .        'Basic'   
---- VAR x3          .             0.5000        +INF             .        'Basic'   

Same primal solution, but different basis. No hint of dual degeneracy however. Bonus question: Which of these two simplex solutions is more correct?

Actually both are correct. From linear programming 101 we know  \(y = c_B^T B^{-1}\) where \(y\) are the duals.  We can calculate this quickly with some R. Suppose we have:

A
##    x0 x1        x2        x3 s1 s2 s3
## e1  1  0 0.3333333 0.6666667  1  0  0
## e2  0  1 0.6666667 0.3333333  0  1  0
## e3  1  1 1.0000000 1.0000000  0  0  1
c
## x0 x1 x2 x3 s1 s2 s3 
##  0  0 -1 -1  0  0  0 

Then we can do:

bs ← c('x2','x3','s2')
t(c[bs]) %*% solve(A[,bs])
##      e1 e2 e3
## [1,]  0  0 -1
bs ← c('x2','x3','s3')
t(c[bs]) %*% solve(A[,bs])
##      e1 e2 e3
## [1,] -1 -1  0

There is even a third optimal basis. It has the following duals

bs ← c('x2','x3','s1')
t(c[bs]) %*% solve(A[,bs])
##      e1 e2 e3
## [1,]  0  0 -1

In a sense we have three optimal solutions and SciPy's simplex solver is not able to find a single one of them.

Conclusion

In general the Simplex solver from scipy.optimize.linprog should be avoided. It fails on quite a few small models. The new interior-point solver available in SciPy 1.0 seems to behave much better.

References

1. Linear programming using scipy.optimize.linprog returns suboptimal solution,  https://stackoverflow.com/questions/46702953/linear-programming-using-scipy-optimize-linprog-returns-suboptimal-solution
2. linprog(method='simplex'), https://docs.scipy.org/doc/scipy/reference/optimize.linprog-simplex.html
3. linprog(method='interior-point'), https://docs.scipy.org/doc/scipy/reference/optimize.linprog-interior-point.html
4. What is this EPS in a GAMS solution?, http://yetanothermathprogrammingconsultant.blogspot.com/2017/09/what-is-this-eps-in-gams-solution.html
5. Finding all optimal LP solutions, http://yetanothermathprogrammingconsultant.blogspot.com/2016/01/finding-all-optimal-lp-solutions.html