Add an interface to SciPy-accessible LP solvers

[michaels0m]

Scipy version 1.6.1 and later versions have access to HiGHS solvers (https://www.maths.ed.ac.uk/hall/HiGHS/#guide) for LPs via the linprog function which are comparable with GLPK and other solvers in terms of performance especially for large sparse LP problems which uses SciPy.sparse matrices.

The aim would be to have a new interface with SciPy so that users can access the HiGHS solvers through CVXPY.

I have forked and cloned the repo and have already started to work on such an interface.

[michaels0m]

I have been working on a new Solver interface to incorporate the HiGHS solvers (https://www.maths.ed.ac.uk/hall/HiGHS/#guide) via SciPy and by following the steps in the contributing section of CVXPY (https://www.cvxpy.org/contributing/index.html#solver-interfaces) I managed to get a working solution however I had two questions:

1. The SciPy linprog() function has a parameter called "method" (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html), however that parameter is already being used in Problem.solve() (line 385 in https://github.com/cvxpy/cvxpy/blob/master/cvxpy/problems/problem.py) . Should I ask users to specify options within a separate dictionary e.g. prob.solve(verbose=True, solver="SCIPY", scipy_options={method='highs-ipm', 'options'={'tol':1e-6}}) ? Or should I change the name of the parameter to something like "scipy_method" and it can then be included in the documentation?

2. linprog() does not return any dual variables, however it does return the slack variables. In the invert function within the scipy_conif.py file I created, I return the dual_vars as None. Are there any consequences of doing this, or is there a way to calculate the dual variables which isn't time consuming?

[rileyjmurray]

For (1), I think having an explicit keyword argument scipy_method is the way to go. That wouldn't go in the docs for the Problem class, it would just be mentioned in this part of the web docs. The SCIPY class should also have the ability to process arguments like

d={'method':'highs-ipm', 'tol':1e-6}
prob.solve(solver='SCIPY', scipy_options=d)

then your interface will effectively do something like

method = d.pop('method', 'interior_point')  # 'interior_point' is the SciPy default
# ^ extract any keyword arguments that SciPy doesn't put in "options"
result = linprog(...,method=method, options=d)

For (2), setting dual_vars=None is fine. As far as I know, there is no simple way to determine optimal dual variables given only optimal primal variables. It's odd that SciPy doesn't return dual variables, given that the underlying algorithms sure do.

[michaels0m]

@rileyjmurray - thank you for your response. I have now amended the scipy_conif script so that it works as you describe above. I did not change the "method" parameter as it does not cause any issues when contained within scipy_options.

It will pop "method" and "bounds" from the dictionary and leave the rest to be used as options. I did not include "callback" or "x0" as neither are supported by the HiGHS solvers and the latter is only applicable for the "revised-simplex" method. I also changed the default solver to "highs" as suggested within the Notes section of the documnetation: https://docs.scipy.org/doc/scipy/reference/optimize.linprog-highs.html#optimize-linprog-highs

I ran a couple of scenarios and it seems to work really well with insignificant overhead from SciPy, it is still a pity that we do not have access to the dual-variables but we could think of a direct interface with HiGHS solvers in the future if needed.

One question I had before moving on to building tests was on the apply() function within the _conif.py files. Since I only implemented SciPy's linprog function, it will only work for LPs just like with GLPK. I took a look at glpk_conif.py (https://github.com/cvxpy/cvxpy/blob/master/cvxpy/reductions/solvers/conic_solvers/glpk_conif.py) and noticed that it calls the cvxopt_confi.py apply() function through super() which has the following code:

constr_map = problem.constr_map inv_data[self.EQ_CONSTR] = constr_map[Zero] inv_data[self.NEQ_CONSTR] = constr_map[NonNeg] + constr_map[SOC] + constr_map[PSD] len_eq = problem.cone_dims.zero

Since we are only looking at LPs, wouldn't that mean that SOC or PSD or ExpCone constraints are not supported? So should the code for the new SciPy interface (and possible a change for GLPK) look more like this:

constr_map = problem.constr_map inv_data[self.EQ_CONSTR] = constr_map[Zero] inv_data[self.NEQ_CONSTR] = constr_map[NonNeg] len_eq = problem.cone_dims.zero

[michaels0m]

Side note: I looked into the scipy wrapper of the HiGHS solvers and it looks like they are planning to carry back the dual-variables from the solver's results, it just is not implemented yet in the final version. This can be seen in the git repo in this file https://github.com/scipy/scipy/blob/master/scipy/optimize/_linprog_highs.py at lines 388-396 and 415-440 at this point in time.

[rileyjmurray]

Responding to a few points, @Michael-git96

@rileyjmurray - thank you for your response. I have now amended the scipy_conif script so that it works as you describe above. I did not change the "method" parameter as it does not cause any issues when contained within scipy_options.
It will pop "method" and "bounds" from the dictionary and leave the rest to be used as options. I did not include "callback" or "x0" as neither are supported by the HiGHS solvers and the latter is only applicable for the "revised-simplex" method. I also changed the default solver to "highs" as suggested within the Notes section of the documentation: https://docs.scipy.org/doc/scipy/reference/optimize.linprog-highs.html#optimize-linprog-highs

All that sounds good.

I ran a couple of scenarios and it seems to work really well with insignificant overhead from SciPy, it is still a pity that we do not have access to the dual-variables but we could think of a direct interface with HiGHS solvers in the future if needed.

Excellent! For dual variables you can create a new GitHub issue after your PR is merged, and that GitHub issue can be linked to an issue on the SciPy repo.

One question I had before moving on to building tests was on the apply() function within the _conif.py files. Since I only implemented SciPy's linprog function, it will only work for LPs just like with GLPK. I took a look at glpk_conif.py (https://github.com/cvxpy/cvxpy/blob/master/cvxpy/reductions/solvers/conic_solvers/glpk_conif.py) and noticed that it calls the cvxopt_confi.py apply() function through super() which has the following code:
constr_map = problem.constr_map inv_data[self.EQ_CONSTR] = constr_map[Zero] inv_data[self.NEQ_CONSTR] = constr_map[NonNeg] + constr_map[SOC] + constr_map[PSD] len_eq = problem.cone_dims.zero
Since we are only looking at LPs, wouldn't that mean that SOC or PSD or ExpCone constraints are not supported?```

You are correct that ExpCone, SOC, PSD, and PowCone3D are all not supported by this solver. Your interface should inherit from ConicSolver. If you look at conic_solver.py you'll find

cvxpy/cvxpy/reductions/solvers/conic_solvers/conic_solver.py

Line 78 in 524df98

SUPPORTED_CONSTRAINTS = [Zero, NonNeg]

which declares the supported constraints. (You'll also find that ConicSolver.format_constraints has code to account for other cones, but that code won't be invoked unless you modify SCIPY.SUPPORTED_CONSTRAINTS.)

[michaels0m]

@rileyjmurray Thank you, I will ensure that only the NonNeg and Zero constraints are used in the apply() function and I checked and indeed SCIPY does inherit from ConicSolver and so SUPPORTED_CONSTRAINTS only include Zero and NonNeg. I can modify scipy_conif and push the changes so they appear in the pull request I submitted.