ENH: linprog in SciPy 1.2+ discussion

[mdhaber]

scipy.optimize.linprog method='simplex' is in much better shape after several bug fixes and adopting the pre/postprocessing routines originally written for interior-point (closed #5400, #6139, #6690, #7044, #7237, #8174, and #8662 via PRs #8909, #8958, #9081, #9096, #9145 - thanks @Kai-Striega!). Essentially all known linprog bugs have been fixed. Consequently, it's time for some enhancements!

To kick things off, I've submitted PR #9263 to add method = 'revised simplex'. The refactoring in #9145 paved the way for all linprog methods to share pre/postprocessing code, making it very easy to work this new method into the framework.

The fun is far from over, however. I opened this issue to record improvement ideas.

• use SuiteSparse in method = 'interior-point'when available
• allow user-specified linear solver (see custom linear solver for linprog #11201) (not important anymore; HiGHS linear solver is Fine)
• add sparse matrix support for method = 'revised simplex' (ENH: optimize: add HiGHS methods to linprog - continued #12043 adds sparse simplex w/ HiGHS)
• for method='revised simplex', improve routine for driving artificial variables out of the basis by calculating only the required entries of the B^-1 A matrix, and update the basis matrix factorization rather than re-calculating from scratch (ENH: optimize: add HiGHS methods to linprog - continued #12043 adds improved simplex)
• add an option to method = 'revised simplex' to perform the (revised) dual simplex method (and potentially include a heuristic for choosing between the two variants) (ENH: optimize: add HiGHS methods to linprog - continued #12043 adds dual simplex)
• implement the Forest-Tomlin update for method='revised simplex' basis matrix factorization (ENH: optimize: add HiGHS methods to linprog - continued #12043 has simplex with more efficient basis matrix updates)
• add a crossover routine (find basic feasible solution given a solution on the optimal face) for method = 'interior-point' (ENH: optimize: add HiGHS methods to linprog - continued #12043 adds interior point with crossover)
• Add a basis-crashing routine for 'method=revised simplex' (solvers added by ENH: optimize: add HiGHS methods to linprog - continued #12043 will have basis crashing evetually)
• include more sophisticated pivoting rules for method = 'revised simplex' (ENH: optimize: add HiGHS methods to linprog - continued #12043 adds simplex with several pivoting rules)
• allow user to provide initial primal- or dual-feasible solution to method = 'revised simplex'
• add meaningful information about the dual of the original problem, especially for sensitivity analysis, to OptimizeResult produced for all methods (see scipy.optimize.linprog(method='simplex') needs to return duals #8836, How to get Lagrange / lambda multipliers out of 'linprog' optimize subroutine in scipy module ? #11848) (added marginals in ENH: linprog HiGHS marginals/sensitivy analysis #13893; still want ranging)
• add meaningful information about the standard form problem to OptimizeResult - along with details about the relationship between the standard form problem and the original problem (I don't think this is important anymore.)
• add a routine for finding ALL basic optimal solutions to a given problem
• add functionality to aid in sensitivity analysis and updates to the linear programming problem (e.g. adding a constraint or changing a cost). For instance, the OptimizeResult object returned by linprog could include a copy of the standard form problem, solution vector, and any other information (basis, dual variables, etc.) to allow for a warm or hot start (ENH: scipy.optimize: object-oriented interface to HiGHS #15915)
• add tool for reading MPS files (I have code with nearly full functionality that could be cleaned up) (File reader and writer for .mps and .lp files (scipy.optimize.linprog) #14593
• add callback support for method = 'interior-point'/'revised simplex'
• improve the presolve routines: (HiGHS presolve is fine) - additional presolve techniques (see WIP: optimize: added (dense) interpolative decomposition redundancy removal to linprog #9783), iterative resolve, sparse pivot-based redundancy removal, heuristics for pivot-based redundancy removal basis crashing, faster algorithm for redundancy removal (this would help), improve SVD method by removing multiple rows at a time when possible
• add automatic scaling routine (see WIP: optimize: add autoscaling to linprog #10495)
• add most Netlib benchmarks (see BENCH: optimize: more comprehensive linprog benchmarking #10762)
• add redundancy removal benchmarks (see BENCH: optimize: more comprehensive linprog benchmarking #10762)
• add support for integer constraints (I already have a fair amount of work done here that has been waiting for a reliable simplex solver, but there is quite a bit more to be done) (Adds Mixed Integer Linear Programming from highs #14455)
• batch linear programming (Batched linear programs #11195)

Please let me know if you are interested in any of these ideas and I'd be happy to share additional information.

[ev-br]

A comment from the peanut gallery: is "revised simplex" a standard term in the literature? Or is it just that it's revised w.r.t. the original scipy simplex code? (In which case, a more descriptive name might be in order?) Also, it's "revised simplex", but "interior-point" with a dash.

[mdhaber]

@ev-br standard term; different style of implementation (of the same algorithm).
Also it's never hyphenated, but 'interior-point' typically is.

[rgommers]

add support for integer constraints (I already have a fair amount of work done here that has been waiting for a reliable simplex solver, but there is quite a bit more to be done)

Constraints only? You were talking before about a MIP solver, which sounds different from this. However if this is not that, then it's unclear to be what "integer constraints" mean without support for integer variables.

[mdhaber]

A MIP is a linear program with some variables constrained to be integers. That's what I meant by integer constraints, but maybe "integrality constraints" would have been more clear.

[rgommers]

Ah okay, got it. I find either term confusing - constraints in this context are normally the s.t. (subject to) statements in optimization DSLs, what we'd call bounds here. Also in Python context I wouldn't call it a constraints, it's simply using integer variables.

[Kai-Striega]

I'd also like to include #9268, which should strengthen the tests for some of the closed issues

[Kai-Striega]

When #9263 is merged could we consider raising the default tolerance for the simplex method? The existing "bugs" were really only issues with the tolerance - not the underlying logic. When discussed before the consensus to keep the tolerance at 1e-12 was based on the simplex method having a (theoretical) tolerance very close to machine accuracy.

Having a revised simplex solver (ideally as the default) will largely make simplex obsolete. Except, maybe, for educational/small problems and seeing every iteration is important.

@mdhaber @rgommers

[rgommers]

Good point. I don't see a problem with loosening the tolerance.

[rgommers]

I've added a note on this in my roadmap update PR:

Many ideas for additional functionality (e.g. integer constraints, sparse
matrix support, performance improvements) in ``linprog``, see gh-9269.

[mdhaber]

Nice to meet you today, @melissawm!

[BinglePan]

@mdhaber @Kai-Striega Can linprog with simplex or revised-simplex method be used to do sensitivity analysis now? It seems that the new version of linprog no longer supports the tableau and basis fields in the callback function. But without tableau and basis it's very hard to do sensitivity analysis whose formulas are based on its entries.

[mdhaber]

Can linprog with simplex or revised-simplex method be used to do sensitivity analysis now?

No : ( Still on the list above.

It seems that the new version of linprog no longer supports the tableau and basis fields in the callback function.

Yes, it looks like this was changed in gh-9145.

But without tableau and basis it's very hard to do sensitivity analysis whose formulas are based on its entries.

I'm curious how you did sensitivity analysis even when you had the tableau. The tableau corresponds with a canonical form linear programming problem generated from the original input problem, which is typically expressed in a more flexible standard form. Did you always provide your problem in a canonical form in hope that the tableau would correspond with your original problem?

When we were refactoring the linprog code to better accommodate the new methods, the relationship between the original standard form input and the canonical form problem used internally got even more complicated (for better solution efficiency). Even if you provided a canonical form problem as input, the presolver might eliminate constraints or variables (according to complicated rules), making the results of the tableau even more difficult to decipher. So when we changed the callback interface, we decided not to add back the tableau information until we could properly document the relationship between the original input problem and the tableau used by simplex.

Adding the tableau back to the OptimizeResult object passed to callback functions is probably one line of code, but we can't see a use for it unless there is some guide to the relationship between the tableau and the original problem. Can you tell us more about your use case? Maybe there is something we're missing.

[BinglePan]

Can linprog with simplex or revised-simplex method be used to do sensitivity analysis now?

No : ( Still on the list above.

It seems that the new version of linprog no longer supports the tableau and basis fields in the callback function.

Yes, it looks like this was changed in gh-9145.

But without tableau and basis it's very hard to do sensitivity analysis whose formulas are based on its entries.

I'm curious how you did sensitivity analysis even when you had the tableau. The tableau corresponds with a canonical form linear programming problem generated from the original input problem, which is typically expressed in a more flexible standard form. Did you always provide your problem in a canonical form in hope that the tableau would correspond with your original problem?

When we were refactoring the linprog code to better accommodate the new methods, the relationship between the original standard form input and the canonical form problem used internally got even more complicated (for better solution efficiency). Even if you provided a canonical form problem as input, the presolver might eliminate constraints or variables (according to complicated rules), making the results of the tableau even more difficult to decipher. So when we changed the callback interface, we decided not to add back the tableau information until we could properly document the relationship between the original input problem and the tableau used by simplex.

Adding the tableau back to the OptimizeResult object passed to callback functions is probably one line of code, but we can't see a use for it unless there is some guide to the relationship between the tableau and the original problem. Can you tell us more about your use case? Maybe there is something we're missing.

I'm doing sensitivity analysis using linprog. In traditional sensitivity analysis, the formulas for cost coefficients and RHS values are based on the entries of the augmented matrix of the canonical form. I've downgraded Scipy to 1.1.0 to get the tableau and basis fields in the callback function. Though I don't know whether it's what I need, I may get some useful information. I did get the tableau. But as you said, I couldn't understand why it looked like this. Maybe I have to simplify the data model and study the core code for simplex.

[mdhaber]

Ah, so even with the tableau you are having trouble, as expected. (After all, that's why we removed the tableau and haven't yet added it back.)

If you're interested in us adding back the tableau, it would help if you could think about what other information you'd need in order to make the tableau useful. How can we communicate to the user which rows and columns of the tableau correspond with which constraints and variables of the original problem? Keep in mind that when converting a problem to canonical form, variables can be added (artificial variables), removed (in presolve when the value is known), or even split into two (unbounded variables); they are also generally shifted to lie in the range [0, oo). Constraints can also be removed (e.g. when they are redundant or can be trivially satisfied) or added (due to bounds on variables). Given all that, if you can help us figure out how to make the tableau for the canonical form useful, it would go a long way toward adding sensitivity analysis features to linprog.

[BinglePan]

Tableau is indeed useful for sensitivity analysis, at least for educational purpose. You know it's time-consuming to do simplex on paper. But linprog could help. You've mentioned many rules for converting an original LP problem to a canonical form. It doesn't matter only if these rules are standard and specified. I think there are always some users who can find the way out.

[mdhaber]

The "rules" are rather complicated; it's more of an algorithm. So I'm afraid that this is impractical to document. @Kai-Striega with that in mind, what are your thoughts on adding tableau into the OptimizeResult object, even if it is very difficult to interpret?

[Kai-Striega]

The "rules" are rather complicated; it's more of an algorithm. So I'm afraid that this is impractical to document.

I agree. Even if they were simpler I don't think they should be documented - if they are improved - we won't have any issues with backwards compatibility.

what are your thoughts on adding tableau into the OptimizeResult object, even if it is very difficult to interpret?

As it stands, the values passed to callback have been transformed back to the users problems while the tableau still represents the internal problem. I can't see being useful. Passing the internal tableau/variables/et al. could be useful as an educational tool - it shows exactly what the simplex solver is doing each iteration. This would be similar to what the solver returned in < 1.2.

I can't tell how useful this will be for sensitivity analysis. Perhaps without presolve the tableau will be similar enough to reverse engineer the remaining changes? @BinglePan if your interested I'll test a few problems and see the tableau's output.

[BinglePan]

@Kai-Striega I've tested linprog in scipy 1.1.0 on a rather simple example (3 non-negative variables and 3 constraints (1 equality and 2 inequalities)). In my mind, the size of the simplex tableau for this LP problem should be 4 (3 constraints + objective function) * 6 ( 3 decision variables + 2 slack/surplus variables + 1 RHS value). But I got a 5 * 9 matrix instead. It seemed that even if the first constraint was an equality, a surplus variable was also added in. In addition, the last constraint looked like it was generated by all the previous constraints. These rules did confuse me. But in the last iteration, the tableau was reduced to 4*6 size as it should be. But I was not sure whether it was the one I wanted due to these rules.

[Kai-Striega]

@BinglePan the 5×9 matrix represents the phase one tableau i.e. it includes the pseudo-objective function (f') and 3 artificial variables. When a solution to the phase one problem was found the extra row and columns were removed.

Between 1.1 and 1.2 the presolve method used by the interior-point solver was added to the simplex solver. This uses more complex rules/transformations/substitutions to transform the user's problem into the problem used internally. This can be switched off by setting presolve=False. When passing the variables to callback the variables are transformed back into the user's problem, however, this is not defined for the tableau.

[rgommers]

@mdhaber time to update this issue after all the Highs improvements?

[mdhaber]

Done. The important remaining items, colidated here for concision:

• add tool for reading MPS files (see File reader and writer for .mps and .lp files (scipy.optimize.linprog) #14593)
• add support for integer constraints (see WIP: HiGHS MIP wrapper mckib2/scipy#28, ENH: optimize: milp: add function for mixed integer linear programming mckib2/scipy#31. SciPy PRs imminent.)
• add more advanced interface for exploring the problem and changes to the problem (see Adds Mixed Integer Linear Programming from highs #14455. If HiGHS makes this sort of thing possible, we'll eventually add it.)

[mckib2]

HiGHS does have warm-start capability (mainly for MIP), but if we can preserve the underlying LP model between runs, in principle I believe we should get that functionality out of the box (see e.g. ERGO-Code/HiGHS#617). Such an interface could be object-oriented or functional, both would require modifying the Cython wrapper to carry the solved/modified LP model for the lifetime of the Python object. I'm not sure what it would look like (I'm sure we could look at examples from other libraries), but it would need to be a separate control path from the existing linprog function to avoid needlessly rewriting all the data that went unchanged run to run. An OO solution might look like this:

from scipy.optimize import MyCoolLPObj
lp = MyCoolLPObj(c, Aub, bub, Aeq, beq, options, ...)  # construct a new LP problem
lp.solve()  # solve for the first time 
lp.c[4] = 4  # make some small changes to the problem
lp.solve()  # solve again, taking advantage of the HiGHS warm-start machinery
sol = lp.getSolution()  # get solution details

A functional interface is also easy to imagine -- just make MyCoolLPObj the input/output parameters to a solve function

---

Another thing to add to wishlist above might be a primal solver. Not sure if there's a tremendous amount of utility (I certainly don't have a use case), but it should be easy to expose HiGHS' primal solver now that it is a little more mature.

[Kai-Striega]

@mdhaber I'm not sure if this is the right place to start the discussion but I'd like to propose adding type hints to the linear programming solvers. If we're looking at the public functions this would only be linprog, linprog_verbose_callback and linprog_terse_callback with the remaining solvers being added at a later date.

I would be happy to have a go at it but I'd need some review as I haven't done much with type hints yet.

[mdhaber]

I have yet to be convinced of their utility, so I wouldn't be able to review : /

[Kai-Striega]

Ok. Unfortunately I'm in a similar position. The reason I wanted to consider adding type hints is that some users really value them and I think it would be professional to supply them with the option.

[mdhaber]

I think it's safe to close this. All of the important items here are tracked in separate issues, and in this case, I don't think we need a meta-issue.