Cvxpy: Linear mixed-integer solver problem is failing with GLPK-MI

[hadware]

Describe the bug
Linear mixed-integer solver problem is failing with GLPK-MI (and used to work in older versions, namely 1.0.25, the one I used to design my problem).

To Reproduce
Here's the smallest test I could come up with. Sorry for the size of the constraint matrices and the "real" values, but I couldn't reproduce that bug with small handcrafted matrices.

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.], 
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.], 
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.], 
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) == 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve()
print(f"Optimal disorder is {optimal_disorder}")

Expected behavior
It should come up with a solution. Adding verbose=True did not print anything useful that might help me solve that bug.

Output

Traceback (most recent call last):
  File "/home/hadware/Code/CoML/test_cvxpy/linsolver_test.py", line 16, in <module>
    optimal_disorder = prob.solve()
  File "/home/hadware/Code/CoML/test_cvxpy/venv/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 396, in solve
    return solve_func(self, *args, **kwargs)
  File "/home/hadware/Code/CoML/test_cvxpy/venv/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 748, in _solve
    self.unpack_results(solution, solving_chain, inverse_data)
  File "/home/hadware/Code/CoML/test_cvxpy/venv/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 1062, in unpack_results
    "Try another solver, or solve with verbose=True for more "
cvxpy.error.SolverError: Solver 'GLPK_MI' failed. Try another solver, or solve with verbose=True for more information.

Version

• OS: Ubuntu 18.04, python 3.6, venv
• CVXPY Version: 1.1.4 (same bug on 1.1.3)
• CVXOPT Version: 1.2.5

[SteveDiamond]

Are you sure it was previously solving with GLPK_MI? We removed the ECOS_BB solver in 1.1, which would have been the default for your problem.

[hadware]

Ah, good point. Indeed, when I ran it on the older version without specifying a solver (and while setting verbose=True), it prints out ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS (among other things).

I guess this means you're right? Which solver should I use for my problem in more recent versions of cvxpy?

[SteveDiamond]

We're pushing SCIP now: https://www.cvxpy.org/install/#install-with-scip-support
It's unfortunate that people now need to install a solver to solve mixed integer problems. The issue is that the solutions from ECOS_BB were often totally wrong.

[hadware]

Ah, darnit. There's no way for GLPK_MI to be fixed to be able to solve this kind of problem? Having to install a solver by hand is indeed a bit troublesome, since I'm building a small evaluation library for annotation tasks, and asking users to install and configure paths for a solver will probably deter a lot of potential users from using it.

EDIT: Sorry, I'm bad and I should feel bad. I didn't read the doc thoroughly enough: if I understood it correctly, to have GLPK working, I need to manually setup its bindings as well?

[rileyjmurray]

If you install CVXOPT, then you automatically have access to GLPK and GLPK_MI. Installing CVXOPT is very easy with pip.

Also: XPRESS community edition is available from pip, and solves mixed integer problems with at most 500 variables.

_{Sent with GitHawk}

[lotuski]

Same problem~ GLPK_MI works for cvxpy-1.0.10, but currently it failed

[lumbric]

Update: see also comment below.

I'm having similar issues with a similar problem, I can also reproduce the problem in the bug description. After a bit of trial and error, I suspect that when using GLPK_MI, equality constraints are not supported or not working correctly. There might be an easy workaround by adding the constraint with inequalities in both directions. At least for your example it seems to spit out the same result as it does with Gurobi. Disclaimer: lots of guess work involved here :)

Environment

>>> cvxpy.__version__
'1.1.3'
>>> cvxopt.__version__

glpk=4.65=he80fd80_1002

$ /opt/gurobi810/linux64/bin/gurobi_cl --version
Gurobi Optimizer version 8.1.0 build v8.1.0rc1 (linux64)

Complete cona environment

Minimal example 1

the original example with making the solver choice explicit via solver=cp.GLPK_MI.

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) == 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)
print(f"Optimal disorder is {optimal_disorder}")

Output:

Traceback (most recent call last):
  File "/home/p/abcvoting/survey/example05-pav.py", line 16, in <module>
    optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 395, in solve
    return solve_func(self, *args, **kwargs)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 747, in _solve
    self.unpack_results(solution, solving_chain, inverse_data)
  File "/opt/miniconda3/envs/abcvoting/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 1055, in unpack_results
    "Try another solver, or solve with verbose=True for more "
cvxpy.error.SolverError: Solver 'GLPK_MI' failed. Try another solver, or solve with verbose=True for more information.

Process finished with exit code 1

Minimal example 2

same as example 1, but replacing the constraint with two inequalities (and adding more output).

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) <= 1, cp.matmul(A, x) >= 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GLPK_MI, verbose=True)

print(f"Optimal disorder is {optimal_disorder}")
print(f"x: {x.value}")

Output:

/opt/miniconda3/envs/abcvoting/bin/python /home/p/abcvoting/survey/example05-pav.py
      0: obj =   0.000000000e+00 inf =   5.000e+00 (5)
      5: obj =   3.360000000e+00 inf =   0.000e+00 (0)
*     9: obj =   3.360000000e+00 inf =   0.000e+00 (0)
+     9: mip =     not found yet >=              -inf        (1; 0)
+     9: >>>>>   3.360000000e+00 >=   3.360000000e+00   0.0% (1; 0)
+     9: mip =   3.360000000e+00 >=     tree is empty   0.0% (0; 1)
Optimal disorder is 3.36
x: [1. 0. 0. 0. 0. 0. 1. 1. 1. 0.]

Process finished with exit code 0

Minimal example 3

Same thing with Gurobi and equality constraint for cross-checking the solution:

import cvxpy as cp
import numpy as np

x = cp.Variable(shape=10, boolean=True)
disorders = np.array([0.36, 0.36, 0.36, 0.36, 1., 1., 1., 1., 1., 1., ])
A = np.array([[1, 0, 0, 0, 0, 1, 0, 0, 0, 0.],
              [0, 1, 0, 0, 0, 0, 1, 0, 0, 0.],
              [0, 0, 1, 0, 0, 0, 0, 1, 0, 0.],
              [0, 0, 0, 1, 0, 0, 0, 0, 1, 0.],
              [1, 1, 1, 1, 1, 0, 0, 0, 0, 0.]])


obj = cp.Minimize(disorders.T @ x)
constraints = [cp.matmul(A, x) == 1]
prob = cp.Problem(obj, constraints)
optimal_disorder = prob.solve(solver=cp.GUROBI, verbose=True)

print(f"Optimal disorder is {optimal_disorder}")
print(f"x: {x.value}")

Output:

Parameter OutputFlag unchanged
   Value: 1  Min: 0  Max: 1  Default: 1
Changed value of parameter QCPDual to 1
   Prev: 0  Min: 0  Max: 1  Default: 0
Optimize a model with 5 rows, 10 columns and 13 nonzeros
Variable types: 0 continuous, 10 integer (10 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [4e-01, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 5.0000000
Presolve removed 5 rows and 10 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 8 available processors)

Solution count 2: 3.36 5 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.360000000000e+00, best bound 3.360000000000e+00, gap 0.0000%
Optimal disorder is 3.36
x: [1. 0. 0. 0. 0. 0. 1. 1. 1. 0.]

Process finished with exit code 0

[hadware]

Woah great work @lumbric ! That's actually quite an hacky yet elegant way to make this work. I'll try to see if my unit tests fail or not when I'm using this workaround.

[SteveDiamond]

@hadware not sure if you know this, but we added ECOS BB back as a solver. You need to call it explicitly (prob.solve(solver=cp.ECOS_BB)).

[hadware]

@SteveDiamond nice! However, in a previous message, you said that

The issue is that the solutions from ECOS_BB were often totally wrong.

Are these issues solved?

[lumbric]

The issue is that the solutions from ECOS_BB were often totally wrong.

Are these issues solved?

Probably not: from the documentation:

there are recurring correctness issues with ECOS_BB

See also changes in 1.1.6 and #1180.

[SteveDiamond]

@hadware yes there are correctness issues, but ECOS BB was working for you in the past, so it might be a good option over GLPK MI.

[lumbric]

The workaround works, but it's actually too much. It seems to me as if there needs to be at least one in-equality constraint, no matter which one. Otherwise cvxopt.glpk.ilp() will raise ValueError: m must be a positive integer. So it seems to be a bug in CVXPY, which occurs only if there are no in-equality constraints. I'll try to prepare a patch. Unfortunately, I couldn't find any documentation of cvxopt.glpk.ilp() except this stackoverflow Q&A. If there is really no documentation, it's not surprising that bugs can occur in corner cases... :shrug: Edit: there is a docstring (identical to the stackoverflow Q&A).

[lumbric]

It's actually pretty simple: cvxopt.glpk.ilp() doesn't support linear programs without inequality constraint, so it's not really a bug. CVXPY could probably pass on the error message to the user in a nicer format, but I wouldn't say this is a bug. But it should be possible to add support for it easily, I'm just not sure what's the right way to go here. Is there any reason, why one would require at least one in-equality constraint for the solver?

I guess we could simply use the equality constraints also for the in-equality constraints if there is no in-equality constraint, i.e. add the workaround from my previous post to CVXPY. Or we could add just one (arbitrary?) in-equality constraint as equality constraint. I guess solving a linear program without any constraints, doesn't make sense, right? So there must be some constraint. Is there a performance or numerical difference which constraint is chosen? Is there any better solution than this? Maybe it needs to be fixed upstream in CVXOPT?

It might be helpful to improve error handling es well. Debugging might be easier, if a an error such as "m must be a positive integer" is passed to the user.

The docstring of cvxopt.glpk.ilp():

Solves a mixed integer linear program using GLPK.

(status, x) = ilp(c, G, h, A, b, I, B)

PURPOSE
Solves the mixed integer linear programming problem

    minimize    c'*x
    subject to  G*x <= h
                A*x = b
                x[k] is integer for k in I
                x[k] is binary for k in B

ARGUMENTS
c            nx1 dense 'd' matrix with n>=1

G            mxn dense or sparse 'd' matrix with m>=1

h            mx1 dense 'd' matrix

A            pxn dense or sparse 'd' matrix with p>=0

b            px1 dense 'd' matrix

I            set of indices of integer variables

B            set of indices of binary variables

I'm new to the internals of CVXPY, but would be happy to provide a patch if somebody could help me a bit along the way.

[SteveDiamond]

Thanks for getting to the bottom of this! You could try adding the constraint 0*x <= 0, but I suspect the solver will fail.

[rileyjmurray]

@SteveDiamond the best solution from CVXPY's standpoint is to tell GLPK our binary variables are merely integer variables, and to add inequalities 0 <= x and x <= 1 to the canonicalized problem data (as appropriate). This could be skipped if the problem already has inequality constraints from the user's side.

[lumbric]

@rileyjmurray oh I think our minimal example wasn't minimal. boolean=True is not necessary to reproduce the problem, so I guess your suggestion wouldn't work in this case:

import cvxpy as cp
import numpy as np

sum = cp.atoms.affine.sum.sum

vars = cp.Variable(5)

constraints = [sum(vars) == 3]

score = sum(cp.atoms.affine.binary_operators.multiply(vars, np.arange(5)))
objective = cp.Maximize(score)

problem = cp.Problem(objective, constraints)

problem.solve(solver=cp.GLPK_MI, verbose=True)

print("-------- Result -----------")

print(problem)
print("\n")
print(vars.value)

[rileyjmurray]

@lumbric yes, but your original example was still a perfectly interesting optimization problem. Users should not call GLPK for problems that are purely linear-algebraic. I think the only interesting case that remains is when you have linear equations and integer constraints but not binary constraints. In that case the only thing to do is formulate the problem with two-sided inequalities. I’ll note that may seem like a bad hack, but back in “the old days” that was a thing people needed to do regularly to get problems into standard form for simplex-based LP solvers.

_{Sent with GitHawk}

[lumbric]

In that case the only thing to do is formulate the problem with two-sided inequalities

Hm yes, that's what I suggested in the first workaround. Then I found out that one of the inequalities would be enough. @SteveDiamond suggested to go for 0*x <= 0.

I wonder if there is any (numerical?) difference, depending on which approach we choose.