Infeasible solution for an lp even though there exists feasible solution(using cvxopt python)

Question

I am trying to find an lp solution to the following problem and even though I can construct feasible points by hand , I seem to get a infeasible certificate from cvxopt. Below is the example and snippet of code

import numpy as np
from cvxopt import solvers, matrix
A = np.array([[0.221, -11.76, -1], [0.2017, -9.88, -1], [0.176, -8.045, 
-1], [0.253, 1.72, 0.]])
print(A)
[[  0.221  -11.76    -1.    ]
[  0.2017  -9.88    -1.    ]
[  0.176   -8.045   -1.    ]
[  0.253    1.72     0.    ]]
b = np.array([-147459573, -139576175, -131383060, 20e6])
x0 = np.array([0.,0.,-np.min(b)])
x1 = np.array([20.e6,8.0e6,71e6])
y0 = np.dot(A,x0)
print(y0<=b)
[ True  True  True  True]
y1 = np.dot(A,x1)
print(y1 <= b)
[ True  True  True  True]
c = np.array([0.,0.,1.])
A = matrix(A)
b = matrix(b)
c = matrix(c)
res = solvers.lp(c,A,b)
pcost       dcost       gap    pres   dres   k/t
0:  1.0031e+08  1.4805e+09  4e+06  2e-03  1e+02  1e+00
Certificate of primal infeasibility found.

As one can see from above x0, x1 clearly are in the feasible set but the solution seems to say that primal is infeasible. Am I looking at this wrong ?

Answer 1

Problem

Your problem is very badly scaled as there are very large and very small coefficients. As all those solvers are working with limited-precision floats, this introduces numerical-instabilities.

The usual approach then is problem scaling or reformulation. There are tons of books and probably papers too (mostly in some chapter about preprocessing), but i'm just citing Mosek's docs here as this is readily available:

Problems containing data with large and/or small coefficients, say 1.0e+9 or 1.0e-7 , are often hard to solve. Significant digits may be truncated in calculations with finite precision, which can result in the optimizer relying on inaccurate calculations. Since computers work in finite precision, extreme coefficients should be avoided. In general, data around the same “order of magnitude” is preferred, and we will refer to a problem, satisfying this loose property, as being well-scaled. If the problem is not well scaled, MOSEK will try to scale (multiply) constraints and variables by suitable constants. MOSEK solves the scaled problem to improve the numerical properties.

The scaling process is transparent, i.e. the solution to the original problem is reported. It is important to be aware that the optimizer terminates when the termination criterion is met on the scaled problem, therefore significant primal or dual infeasibilities may occur after unscaling for badly scaled problems. The best solution to this problem is to reformulate it, making it better scaled.

By default MOSEK heuristically chooses a suitable scaling. The scaling for interior-point and simplex optimizers can be controlled with the parameters MSK_IPAR_INTPNT_SCALING and MSK_IPAR_SIM_SCALING respectively.

Painless alternative (at least for this problem)

Although ecos (conic solver; open-source) is ready to solve much more complex problems, it seems to do much better preprocessing here and can solve your problem. The modelling-framework which is calling ecos is cvxpy:

import numpy as np

A = np.array([[0.221, -11.76, -1],
              [0.2017, -9.88, -1],
              [0.176, -8.045, -1],
              [0.253, 1.72, 0.]])

b = np.array([-147459573, -139576175, -131383060, 20e6])
c = np.array([0, 0, 1])

""" TRY CVXPY + ECOS """

from cvxpy import *
x = Variable(3)
constraints = [A*x <= b, x >= 0]
objective = Minimize(c * x)
problem = Problem(objective, constraints)
problem.solve(verbose=True)
print('problem state: ', problem.status)
print('solution: ')
x_sol = np.array(x.value.flat)
print(x_sol)
print( (A.dot(x_sol) <= b) )

Output:

ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +7.905e+06  +1.394e+08  +1e+08  2e-01  9e-01  1e+00  2e+07    ---    ---    1  1  - |  -  -
 1  +1.415e+07  +5.847e+07  +3e+07  7e-02  3e-01  6e+05  4e+06  0.7593  3e-02   0  0  0 |  0  0
 2  +3.447e+07  +5.278e+07  +1e+07  3e-02  2e-01  1e+06  2e+06  0.8506  2e-01   0  0  0 |  0  0
 3  +4.015e+07  +4.038e+07  +1e+06  4e-04  9e-02  2e+04  2e+05  0.9890  1e-03   0  0  0 |  0  0
 4  +3.820e+07  +3.820e+07  +1e+05  5e-06  4e-02  3e+02  2e+04  0.9856  1e-04   0  0  0 |  0  0
 5  +3.784e+07  +3.784e+07  +3e+03  6e-08  8e-03  7e+01  4e+02  0.9890  2e-03   0  0  0 |  0  0
 6  +3.784e+07  +3.784e+07  +4e+01  7e-10  2e-04  9e-01  5e+00  0.9890  1e-04   0  0  0 |  0  0
 7  +3.784e+07  +3.784e+07  +4e-01  8e-12  5e-06  1e-02  5e-02  0.9890  1e-04   0  0  0 |  0  0
 8  +3.784e+07  +3.784e+07  +4e-03  9e-14  7e-08  1e-04  6e-04  0.9890  1e-04   0  0  0 |  0  0
 9  +3.784e+07  +3.784e+07  +5e-05  1e-15  1e-09  2e-06  6e-06  0.9890  1e-04   1  0  0 |  0  0

OPTIMAL (within feastol=1.0e-09, reltol=1.3e-12, abstol=5.0e-05).
Runtime: 0.003211 seconds.

problem state:  optimal
solution:
[  5.49533698e-05   1.16279070e+07   3.78365484e+07]
[ True  True  True  True]

Mosek should be able to solve this too!