MixedIntegerLinearProgram/HybridBackend: Reconstruct exact rational/algebraic basic solution

[mkoeppe]

Sometimes one can use a fast numerical LP solver to solve a problem to "optimality", then reconstruct the primal and dual solution in rational arithmetic (or over whatever base_ring was used...) and in this way prove that this basis is indeed optimal. MixedIntegerLinearProgram should support this mode of operation.

The current branch, on top of #20296, attempts to do this by implementing a HybridBackend, which delegates to two backends:

• a fast, possibly inexact backend (Gurobi or GLPK or even GLPK with glp_exact -- see #18764)
• a slow, exact one that can set the simplex basis (only InteractiveLPBackend fits the bill - from #20296)

Ideally, in pure LP mode, both backends would support the basis-status functions that can transplant the (hopefully) optimal (hopefully-)basis from the inexact LP to the exact LP.

If the inexact LP cannot provide a basis (because its "basis" is not a basis due to numerics, or because basis-status functions are not available), one could at least try to make use of the numerical solution vector and try to reconstruct a basis, like in interior-point-to-simplex crossover (a classical paper: http://www.caam.rice.edu/caam/trs/91/TR91-32.pdf)

In MIP mode, could at least try to set the cleaned-up numerical solution vector as a known solution, to speed up branch-and-cut in the exact solver.

Sounds like a big ticket; we'll do this step by step.

18685 provides the necessary basis-status functions (for the GLPK backend).

18688 provides a solver-independent interface to these functions.

18804 exposes basis status via backend dictionaries.

Depends on #18685 Depends on #18688 Depends on #20296

CC: @yuan-zhou @nathanncohen @dimpase

Component: numerical

Author: Matthias Koeppe, Yuan Zhou

Branch/Commit: u/yzh/hybrid_backend @ 50773ff

Issue created by migration from https://trac.sagemath.org/ticket/18735

[dimpase]

comment:2

Is ppl (pplLP backend, which works with exact arithmetic) too slow for you?

[dimpase]

comment:3

On the other hand, a solver-independent way to get an optimal dual solution is very much welcome, as this is lacking currently, and often needed.

[dimpase]

Description changed:

--- 
+++ 
@@ -1,6 +1,6 @@
 Sometimes one can use a fast numerical LP solver to solve a problem to "optimality", 
 then reconstruct the primal and dual solution in rational arithmetic (or over whatever base_ring was used...) and in this way prove that this basis is indeed optimal. 
-MixedIntegerLinearProgram should support this mode of operation.
+`MixedIntegerLinearProgram` should support this mode of operation.

 #18685 provides the necessary basis-status functions (for the GLPK backend).
 #18688 provides a solver-independent interface to these functions.

[mkoeppe]

Description changed:

--- 
+++ 
@@ -2,6 +2,8 @@
 then reconstruct the primal and dual solution in rational arithmetic (or over whatever base_ring was used...) and in this way prove that this basis is indeed optimal. 
 `MixedIntegerLinearProgram` should support this mode of operation.

+This would be particularly interesting in conjunction with #18764. (But see #18765 for a different approach.)
+
 #18685 provides the necessary basis-status functions (for the GLPK backend).
 #18688 provides a solver-independent interface to these functions.

[mkoeppe]

comment:5

Replying to @dimpase:

Is ppl (pplLP backend, which works with exact arithmetic) too slow for you?

Dima, ppl's implementation of the double description method is very good, but its LP solver is not suitable for problems of even moderate sizes.

[dimpase]

comment:6

Replying to @mkoeppe:

Replying to @dimpase:

Is ppl (pplLP backend, which works with exact arithmetic) too slow for you?

Dima, ppl's implementation of the double description method is very good, but its LP solver is not suitable for problems of even moderate sizes.

Would you mind providing an example of PPL choking on an LP doable in exact arithmetic by another solver? We use PPL's LP solver in codesize_upper_bound(...,algorithm="LP") and never saw a problem... (Although perhaps the difficulty from entry sizes dominate the the one from the dimension in this case).

[mkoeppe]

comment:7

Replying to @dimpase:

Would you mind providing an example of PPL choking on an LP doable in exact arithmetic by another solver? We use PPL's LP solver in codesize_upper_bound(...,algorithm="LP") and never saw a problem... (Although perhaps the difficulty from entry sizes dominate the the one from the dimension in this case).

In our experiments here, we don't actually have numerical difficulties with floating-point based solvers; we just want to be sure that we have an exact optimal solution. With #18764 (glp_exact; please review) we have now run some tests to compare performance:

                                glp_simplex                glp_simplex+glp_exact
   glp_simplex    glp_exact     +glp_exact    ppl          + reconstruction in Sage
10  4.20            51.92             7.78    207.07          289.00
11  5.08            58.49             9.43    3451.42         574.72
12  7.55           101.72            11.32    1252.91         808.73
13  7.21           279.08            13.57    1424.28        1019.95
14  8.41           562.97            15.91    7343.37        1628.54
15 13.10           550.46            18.48    3667.93        2550.94

As you can see, PPL is much slower than pure glp_exact, and orders of magnitudes slower than glp_simplex followed by glp_exact.

However, currently when we try to reconstruct the solution from the combinatorial basis information, Sage's super slow matrix functions over the rationals get us back to roughly the same order of magnitude as PPL.

It would be interesting to know how the solvers perform on the kind of LPs that you have in mind.

[dimpase]

comment:8

Replying to @mkoeppe:

It would be interesting to know how the solvers perform on the kind of LPs that you have in mind.

LPs I get would be not possible to even enter into a solver without long integers/rationals. That's e.g. behind this function call:

sage:  codesize_upper_bound(70,8,2,algorithm="LP")
9695943911863423

more explicitly, you can do

sage: v,p,r=delsarte_bound_hamming_space(70,8,2,return_data=True)
sage: p
Mixed Integer Program  ( maximization, 71 variables, 148 constraints )

constrains of p have entries as big as 112186277816662845432.

[mkoeppe]

Description changed:

--- 
+++ 
@@ -7,3 +7,5 @@
 #18685 provides the necessary basis-status functions (for the GLPK backend).
 #18688 provides a solver-independent interface to these functions.

+The reconstructed solution could be presented via #20296.
+

[mkoeppe]

Description changed:

--- 
+++ 
@@ -2,10 +2,11 @@
 then reconstruct the primal and dual solution in rational arithmetic (or over whatever base_ring was used...) and in this way prove that this basis is indeed optimal. 
 `MixedIntegerLinearProgram` should support this mode of operation.

-This would be particularly interesting in conjunction with #18764. (But see #18765 for a different approach.)
+The current branch, on top of #20926, attempts to do this by implementing a `HybridBackend`, which delegates to two backends:
+- a fast, possibly inexact backend (Gurobi or GLPK or even GLPK with glp_exact -- see #18764) 
+- a slow, exact one that can set the simplex basis (only `InteractiveLPBackend` fits the bill - from #20296)

 #18685 provides the necessary basis-status functions (for the GLPK backend).
 #18688 provides a solver-independent interface to these functions.
+#18804 exposes basis status via backend dictionaries.

-The reconstructed solution could be presented via #20296.
-

[mkoeppe]

Changed dependencies from #18685, #18688 to #18685, #18688, #20296

[mkoeppe]

Description changed:

--- 
+++ 
@@ -2,7 +2,7 @@
 then reconstruct the primal and dual solution in rational arithmetic (or over whatever base_ring was used...) and in this way prove that this basis is indeed optimal. 
 `MixedIntegerLinearProgram` should support this mode of operation.

-The current branch, on top of #20926, attempts to do this by implementing a `HybridBackend`, which delegates to two backends:
+The current branch, on top of #20296, attempts to do this by implementing a `HybridBackend`, which delegates to two backends:
 - a fast, possibly inexact backend (Gurobi or GLPK or even GLPK with glp_exact -- see #18764) 
 - a slow, exact one that can set the simplex basis (only `InteractiveLPBackend` fits the bill - from #20296)

[mkoeppe]

Branch: u/mkoeppe/hybrid_backend

[mkoeppe]

Commit: 0b8b78a

[mkoeppe]

Description changed:

--- 
+++ 
@@ -6,6 +6,14 @@
 - a fast, possibly inexact backend (Gurobi or GLPK or even GLPK with glp_exact -- see #18764) 
 - a slow, exact one that can set the simplex basis (only `InteractiveLPBackend` fits the bill - from #20296)

+Ideally, in pure LP mode, both backends would support the basis-status functions that can transplant the (hopefully) optimal (hopefully-)basis from the inexact LP to the exact LP.
+
+If the inexact LP cannot provide a basis (because its "basis" is not a basis due to numerics, or because basis-status functions are not available), one could at least try to make use of the numerical solution vector and try to reconstruct a basis, like in interior-point-to-simplex crossover (a classical paper: http://www.caam.rice.edu/caam/trs/91/TR91-32.pdf)
+
+In MIP mode, could at least try to set the cleaned-up numerical solution vector as a known solution, to speed up branch-and-cut in the exact solver. 
+
+Sounds like a big ticket; we'll do this step by step.
+
 #18685 provides the necessary basis-status functions (for the GLPK backend).
 #18688 provides a solver-independent interface to these functions.
 #18804 exposes basis status via backend dictionaries.

[mkoeppe]

Last 10 new commits:

e2319b5	InteractiveLPBackend.get_variable_value: Guard against standard-form transformations
e27f297	InteractiveLPBackend: Make base_ring an init argument
5b0954f	InteractiveLPBackend._variable_type_from_bounds: Add doctests
c4b93aa	InteractiveLPBackend: Fix old-style raise statements
b0a3c1c	GenericBackend: Add a missing '# optional - Nonexistent_LP_solver'
3770be0	default_mip_solver: Handle 'InteractiveLP'
d91c776	default_mip_solver, get_solver: Mention InteractiveLP in the documentation
eaede28	get_solver: Add optional base_ring argument
184249d	MixedIntegerLinearProgram: New base_ring init argument
0b8b78a	HybridBackend: first draft

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0b8b78a to 26dad94

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

26dad94

HybridBackend: first draft

[yuan-zhou]

Changed branch from u/mkoeppe/hybrid_backend to u/yzh/hybrid_backend

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 26dad94 to 5ee7738

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

5ee7738

change Cython to Python code

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

3431fc4

HybridBackend heritage of GenericBackend

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 5ee7738 to 3431fc4

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 3431fc4 to c1951c9

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

c1951c9

fix docstring style

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

e176757	add method objective_constant_term()
c53b4e0	make iterable coefficients a list so that each backend can loop over coefficients
c15aaac	fix doctests

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c1951c9 to c15aaac

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c15aaac to 44bcb39

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

3e4e7e9	fix docstrings
69a9775	warm-start interactive_backend.solve() by providing basic_variables
676c74f	warm-start LPProblemStandardForm.run_revised_simplex_method by providing basic_variables
44bcb39	HybridBackend.solve() reconstruct exact solution

[yuan-zhou]

Author: Yuan Zhou

[yuan-zhou]

Changed author from Yuan Zhou to Matthias Koeppe, Yuan Zhou

[mkoeppe]

comment:25

I don't like the changes to src/sage/numerical/interactive_simplex_method.py too much.

Instead of extending InteractiveLPProblemStandardForm.run_revised_simplex_method, I think it's better to use the "run_simplex_method" method of the current_dictionary (a revised dictionary) directly.

Same amount of code, but would avoid making interactive_simplex_method more complicated.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

5b19c41	Revert "warm-start LPProblemStandardForm.run_revised_simplex_method by providing basic_variables"
1f9d00a	revise interactive_backend.solve() given starting basis

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 44bcb39 to 1f9d00a

[yuan-zhou]

comment:27

Replying to @mkoeppe:

I don't like the changes to src/sage/numerical/interactive_simplex_method.py too much.

Instead of extending InteractiveLPProblemStandardForm.run_revised_simplex_method, I think it's better to use the "run_simplex_method" method of the current_dictionary (a revised dictionary) directly.

Same amount of code, but would avoid making interactive_simplex_method more complicated.

Good idea. Done.

[mkoeppe]

comment:28

Also, I think it would be better to avoid changing the solve interface:

-    cpdef int solve(self) except -1:
+    cpdef int solve(self,  basic_variables=[]) except -1:

Setting an initial basis could be done in a new method perhaps called warmstart(basic_variables) or set_basis(basic_variables)

[mkoeppe]

comment:29

And note that generic backend already defines getter methods for the basis status is_variable_basic, is_variable_nonbasic_at_lower_bound, is_slack_variable_basic, is_slack_variable_nonbasic_at_lower_bound - all of which use variable/row indices.

So this interface:

+            ``basic_variables`` can be one of the following:
+
+            - a list of indices. The indices (starting at 1) correspond to that of the vector formed by `self.interactive_lp_problem().decision_variables()` and `self.interactive_lp_problem().slack_variables()`. Remark that `self.interactive_lp_problem()` can have more variables and constraints than that of `self` if `self` has free variables or `==` constraints.
+
+            - a list of the names of the variables in `self.interactive_lp_problem()`.

looks a bit out of place.

[yuan-zhou]

comment:30

What would be a better input form? Remove the second option and take only a list of indices?

It is tricky because the indices here correspond to the variables in the interactive_lp_problem().standard_form(), and they not necessarily the same as the indices for is_variable_basic etc. which correspond to interactive_lp_problem(). See the following example.

from sage.numerical.backends.generic_backend import get_solver
h = get_solver(solver = ("InteractiveLP"))
h.add_variables(1, lower_bound=None, upper_bound=None);
h.add_variables(1, lower_bound=0, upper_bound=None);
h.add_linear_constraint([(0,2),(1,-1)],-1,None)
h.add_linear_constraint([(0,1),(1,-1)],None, 1)
h.add_linear_constraint([(0,1),(1,1)],2,2)
h.set_objective([0,-1])
lp = h.interactive_lp_problem(); view(lp)
st = lp.standard_form(); view(st)
R = st.coordinate_ring(); R

Replying to @mkoeppe:

And note that generic backend already defines getter methods for the basis status is_variable_basic, is_variable_nonbasic_at_lower_bound, is_slack_variable_basic, is_slack_variable_nonbasic_at_lower_bound - all of which use variable/row indices.

So this interface:

+            ``basic_variables`` can be one of the following:
+
+            - a list of indices. The indices (starting at 1) correspond to that of the vector formed by `self.interactive_lp_problem().decision_variables()` and `self.interactive_lp_problem().slack_variables()`. Remark that `self.interactive_lp_problem()` can have more variables and constraints than that of `self` if `self` has free variables or `==` constraints.
+
+            - a list of the names of the variables in `self.interactive_lp_problem()`.

looks a bit out of place.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

dc36f82

set InteractiveLPBackend.current_dictionary and warm start solve from there

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 1f9d00a to dc36f82

[yuan-zhou]

comment:33

Replying to @mkoeppe:

Also, I think it would be better to avoid changing the solve interface:

-    cpdef int solve(self) except -1:
+    cpdef int solve(self,  basic_variables=[]) except -1:

Setting an initial basis could be done in a new method perhaps called warmstart(basic_variables) or set_basis(basic_variables)

InteractiveLPBackend now has a new method set_dictionary and a new attribute current_dictionary. I didn't name it set_basis, to avoid confusion between the basic variables in a dictionary and those of the generic backend.

[mkoeppe]

comment:34

I still have concerns about the new interface as well. The problem now is that set_basis is not defined for any backend other than InteractiveLPBackend (and also its interface refers to internal details of that).

But more importantly, there should be at least one example that illustrates that this new backend does something useful.

[yuan-zhou]

comment:35

Replying to @mkoeppe:

I still have concerns about the new interface as well. The problem now is that set_basis is not defined for any backend other than InteractiveLPBackend (and also its interface refers to internal details of that).

Do you mean set_dictionary? It has a different interface than set_basis (in #18688 to be written) since the interactive_simplex_method is indeed very special in running the simplex method on the standard form where more auxiliary variables are introduced. I couldn't think of a way to unify InteractiveLPBackend.set_dictionary with other backends' set basis methods.

But more importantly, there should be at least one example that illustrates that this new backend does something useful.

True. The current implementation enables what the first paragraph of this ticket states, by using solver = ("GLPK", "InteractiveLP") as in the examples in solve. However, it is hard to argue with doctests that this new hybrid backend is useful. Shall we close the ticket with the tag "wontfix" as in #18804?

[mkoeppe]

comment:36

Replying to @yuan-zhou:

Replying to @mkoeppe:

I still have concerns about the new interface as well. The problem now is that set_basis is not defined for any backend other than InteractiveLPBackend (and also its interface refers to internal details of that).

Do you mean set_dictionary?

Yes, that's what I meant.

[mkoeppe]

comment:37

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dc36f82 to 71ae94d

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

71ae94d

Replace not hasattr(self, 'xyz') by self.xyz is None

[mkoeppe]

comment:39

Setting a new milestone for this ticket based on a cursory review.

[dimpase]

comment:40

please rebase

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 71ae94d to b6dd9ab