add option to change strategy for solving first root LP in MIP solver

[lgottwald]

Simple change that introduces more options to solve the first LP in the root node of a MIP. Added options to use IPM as well as parallel simplex.

[jajhall]

Thanks. I'll have a look at this, and add a unit test

[lgottwald]

I think there is one issue left. There is currently a simplex iteration limit set based on number of iterations used for the first LP. Need to see that the number of iterations for the first LP is either set to an estimated equivalent of the amount of work used by IPM, or that the limit is not set for the second LP if IPM was used for the first one.

[lgottwald]

Otherwise the lit might be set too low and when the iteration limit is reached IPM will be used again. Once an LP is solved to optimality after the root LP then the average iterations of the LPs solved so far will be used with some factor for the limit, so that it can settle in from then.

[lgottwald]

Also the limit is never set below 10000 so the issue will only appear for harder LPs.

[jajhall]

I know that we'll need it if we are to have deterministic concurrent solution of the root node, but creating synthetic estimates of the work required to solve an LP by simplex and IPM is somewhat non-trivial, even if the important decisions can be made with very inaccurate measures.

I'm not sure exactly what you're suggesting, as I'd have thought that the number of iterations required to solve the root node - without any advanced basis - is likely to be very much more than the number required to solve an LP when there is even an alien basis. That said, I know that you know this, and you're in a better position than I am to have a feel for the actual difference.

[lgottwald]

Yes, but what I mean here is that this change already has some side-effects because the number of simplex iterations for the first root LP will be zero and this number is used to determine iteration limits for later LP solves. It can as a consequence now happen that the limits are set to too strict values and no LP can ever be solved with simplex for the remaining MIP solve. I think I will add a simple workaround to not use the iteration limit when no simplex solves have been made so far.

[lgottwald]

There is some error when solving the first LP without simplex and then switching to simplex. Not everything seems to be initialized in all cases:

==34843== Conditional jump or move depends on uninitialised value(s)
==34843==    at 0x4AE452E: HighsSimplexAnalysis::simplexTimerStart(int, int) (src/simplex/HighsSimplexAnalysis.cpp:596)
==34843==    by 0x4AAE26A: HEkk::computeFactor() (src/simplex/HEkk.cpp:2070)
==34843==    by 0x4AADF55: HEkk::initialiseSimplexLpBasisAndFactor(bool) (src/simplex/HEkk.cpp:1561)
==34843==    by 0x4940BF8: formSimplexLpBasisAndFactor(HighsLpSolverObject&, bool) (src/lp_data/HighsSolution.cpp:1069)
==34843==    by 0x4920E4F: Highs::getBasicVariablesInterface(int*) (src/lp_data/HighsInterface.cpp:1050)
==34843==    by 0x498129D: HighsTableauSeparator::separateLpSolution(HighsLpRelaxation&, HighsLpAggregator&, HighsTransformedLp&, HighsCutPool&) (src/mip/HighsTableauSeparator.cpp:55)
==34843==    by 0x498110F: HighsSeparator::run(HighsLpRelaxation&, HighsLpAggregator&, HighsTransformedLp&, HighsCutPool&) (src/mip/HighsSeparator.cpp:34)
==34843==    by 0x498041F: HighsSeparation::separationRound(HighsDomain&, HighsLpRelaxation::Status&) (src/mip/HighsSeparation.cpp:114)
==34843==    by 0x4954CEF: HighsMipSolverData::rootSeparationRound(HighsSeparation&, int&, HighsLpRelaxation::Status&) (src/mip/HighsMipSolverData.cpp:1086)
==34843==    by 0x49564A9: HighsMipSolverData::evaluateRootNode() (src/mip/HighsMipSolverData.cpp:1350)
==34843==    by 0x494930C: HighsMipSolver::run() (src/mip/HighsMipSolver.cpp:122)
==34843==    by 0x49F121A: HighsPrimalHeuristics::solveSubMip(HighsLp const&, HighsBasis const&, double, std::vector<double, std::allocator<double> >, std::vector<double, std::allocator<double> >, int, int, int) (src/mip/HighsPrimalHeuristics.cpp:124)

This leads to a segfault on roll3000

[lgottwald]

Options for roll3000 segfault are mip_root_lp_strategy = 3 and mip_rel_gap = 0.0

[jajhall]

I can't reproduce this. Runs OK with/without HIGHSINT64

[jajhall]

Behaviour of HiGHS is very odd. I was wanting to explore the execution flow to see if I could understand why you might have observed what you did, so I put a breakpoint at 1232 in HighsMipSolverData.cpp, and one at the start of both the IPX and the simplex solver. I stepped from 1232, and suddenly a different thread stopped at the start of IPX - because it was doing an analytic centre calculation. I've no idea how to inspect the running of the MIP solver

[lgottwald]

You can set the number of threads to 1, then no jumping occurs. In the debugger you can also switch between threads. I will have a look if I can simplify steps to reproduce and can find out more about what happens on my computer.

[lgottwald]

I had a look again and could reproduce the segfault with options either threads=4, threads=2. Different number of threads did not always cause the segfault. But it is also clear from the trace of valgrind I posted what happens. The field analyse_simplex_time in HighsSImplexAnalysis is not initialized. It seems to only be set when setupSimplexTime is called. In the trace I see that the INVERT in this case is triggered by calling getBasicVariablesInterface(), which is weird, but apparently then the uninitialized timer array gets accessed. I call the getBasicVariables in src/mip/HighsTableauSeparator.cpp:55, which should only happen when the LP was solved. I would have expected that calling this function never triggers an invert.

[jajhall]

If getBasicVariables is called and there is no INVERT, then one will be computed for the current basis. I'll try setting threads=2 and putting a breakpoint in src/mip/HighsTableauSeparator.cpp:55 to see if I can understand why no INVERT exists.

[lgottwald]

Could it be that after solving the root LP with IPM no invert may exist?

[jajhall]

Only if crossover isn't performed

[lgottwald]

So I checked this again and on roll3000 still the following happens directly for the very first LP solved with mip_root_lp_strategy = 3:

1. The solver is set to Ipm for the first LP and it is solved
2. In IpxWrapper.cpp Highs goes into the if (have_basic_solution) condition with have_basic_solution = true and imprecise_solution = false
3. ipxBasicSolutionToHighsBasicSolution sets the solution and basis, but leaves basis.alien = true, as it was initialized while also setting basis.valid = true
4. solveLpIpx sets the model status to optimal with valid basis and return kOk
5. In solveLp in HighsSolve.cpp:73 the if condition after running Ipx evaluates to false, solveLp returns without and invertible representation and a valid basis with alien = true

So the first bug is that Highs terminates with kBasisValidityValid while there is no invertible representation. With a valid basis indicated in the info it is expected by the MIP solver that the returned solution is obtained from a valid basis for which Highs has the invertible representation. The second bug is that alien = true after ipxBasicSolutionToHighsBasicSolution.

Following this first two bugs, more bugs are revealed in other functions. When the tableau separator is called it tries to retrieve a tableau row. Highs detects that the invertible representation is missing in getBasicVariablesInterface() and calls formSimplexLpBasisAndFactor().

Now what happens differs on whether alien is true or not, i.e. whether the second bug is fixed. If the bug is not fixed and alien is true, then there is an assertion failure in formSimplexLpBasisAndFactor() because of the value "only_from_known_basis". If the assertion failure is skipped, e.g. in optimized mode or by commenting it out, then formSimplexLpBasisAndFactor() proceeds to go into the if condition at HighsSolution.cpp:1049 to call accommodateAlienBasis, which is the third bug: formSimplexLpBasisAndFactor() returns kOk but ekk_instance.has_invert=false still holds. The "alien" basis passes through the accommodateAlienBasis function fine and needs no accomodation as there is no rank deficiency, but the factor is not retained and is lost once the accommodateAlienBasis function returns. Following the call to accommodateAlienBasis(), formSimplexLpBasisAndFactor() just returns kOk even though there is no factor. I think accommodateAlienBasis should probably check whether the obtained factor has rank_deficiency = 0 and in that case retain it within the ekk_instance of the solver object. In either case formSimplexLpBasisAndFactor() should never return kOk if ekk_instance.has_invert = false still holds. This leads to a segfault in following code.

If the second bug is fixed and alien = false is set inside ipxBasicSolutionToHighsBasicSolution, then there is a different problem. Since simplex has never been called the simplex data is not initialized yet. When HiGHS goes into the same path as in the previous paragraph now formSimplexLpBasisAndFactor() does not call accommodateAlienBasis() and return, but actually does the calls to form the invertible representation by calling initialiseSimplexLpBasisAndFactor(). This uses uninitialized data now, because usually when simplex has been called before, then HEkk::solve() will execute

  initialiseAnalysis();
  initialiseControl();

in lines HEkk.cpp:1033 to initialize clocks before calling initialiseSimplexLpBasisAndFactor() via initialiseForSolve. Now the fourth bug is that initialiseSimplexLpBasisAndFactor() is called directly and uses the uninitialized analysis_ object in computeFactor() to call analysis_.simplexTimerStart(InvertClock).

Of course the other problems disappear when I change the if condition in solveLp in HighsSolve.cpp:73 to always call simplex cleanup when solver_object.basis_.valid_ = true, but I leave it to you to decide whether there is a more elegant solution.

[jajhall]

The setting basis.valid = true has become confusing now that there is the alien distinction. A valid basis is one that has the right number of basic and nonbasic variables.

When IPM+crossover is used to solve a one-off LP, it's more than reasonable to say that the final basis is valid, since since crossover will have performed simplex-like iterations and (I think) the primal and dual values are computed as in the simplex method, rather than being update artifacts of the IPM solver.

However, when solving the root node of a MIP, you are (reasonably) expecting a basic solution as if the LP had been solved by simplex. In particular, you want a basis that is not alien and has an associated invertible representation. I could achieve this by calling the simplex solver to "clean up" after IPM+crossover, but this will have a performance hit for numerically difficult problems where there are primal/dual infeasibilities after factoring the basis and meaningful numbers of simplex iterations are required.

A solution is to call run() a second time after setting solver = "simplex". Normally there will be no iterations. You'll see where I've made this change in HighsLpRelaxation.cpp. Unsurprisingly, the problem with roll3000 disappears.

Clearly this doesn't avoid the assert if someone uses IPM+crossover and then calls Highs::getBasicVariables. It's possible to allow accommodateAlienBasis to be used when only_from_known_basis is true, and return an error if the basis is rank deficient.

An alternative is to perform the simplex clean-up whenever IPM+crossover is used. It has an appeal in that HiGHS would be in the same state when it's found a basic solution to an LP, but I fear for the occasional expensive overhead. I can put it in as an option and do some experiments to decide whether it should always be done.

[lgottwald]

Would it be possible to simply call simplex with a zero iteration limit or to just create the invertible representation without further iterations?

I wonder how it is possible to obtain a basic solution including the values from crossover when no invertible representation is available. Does crossover somehow maintain this internally separate from simplex?

I think when crossover is enabled, then there should always be an invertible representation of the basis, as the user requested to have a basic solution. If that is too expensive for default behavior then I would suggest suggest to set the default to "crossover = false" so that users do not automatically pay the crossover and cleanup overhead.

[lgottwald]

And particularly, I expect that when info says kBasisStatusValid, then the primal and dual solution are the ones that are obtained from inverting the basis given by getBasis() and solving for the values and that this is a valid invertible basis. What is the point otherwise in requesting a basic solution with crossover = true? If I don't need this then I would not even want to do any crossover at all.

[jajhall]

Would it be possible to simply call simplex with a zero iteration limit or to just create the invertible representation without further iterations?

That would create an invertible representation, but either the primal and dual solutions wouldn't correspond to it, or there's a chance that there would be primal/dual infeasibilities after IPX+crossover has determined optimality.

I wonder how it is possible to obtain a basic solution including the values from crossover when no invertible representation is available. Does crossover somehow maintain this internally separate from simplex?

Crossover has its own invertible representation. That's what basiclu is all about. However, after crossover is complete, IPX does things like undo its own scaling - doing some correction - so there's no residual invertible representation that can be inherited.

I'm experimenting with replacing basiclu by HFactor in IPX since I think that I can improve performance. If so, and as I get a better understanding of IPX, it may be possible for the crossover iterations to be performed as if by the simplex solver. Then, when a basis is obtained it's just like a simplex basis. However, IPX is a very sophisticated code, so I'm wary of making modifications without understanding it well.

I think when crossover is enabled, then there should always be an invertible representation of the basis, as the user requested to have a basic solution. If that is too expensive for default behavior then I would suggest suggest to set the default to "crossover = false" so that users do not automatically pay the crossover and cleanup overhead.

I agree, and maybe the overhead of simplex clean-up is something that we just have to accept if we're to be compared with the other IPM solvers that have crossover. And that overhead is likely to diminish if IPX switches to using HFactor rather than basiclu, and as I get a better understanding of IPX.

Also, now that a dual solution is returned when crossover isn't performed, and postsolve without a basis also yields primal and dual solutions, the users who don't need a basic solution can run IPX without crossover. I wouldn't make it the default, though.

And particularly, I expect that when info says kBasisStatusValid, then the primal and dual solution are the ones that are obtained from inverting the basis given by getBasis() and solving for the values and that this is a valid invertible basis. What is the point otherwise in requesting a basic solution with crossover = true? If I don't need this then I would not even want to do any crossover at all.

The point is that, at present, if crossover is run then it yields a vertex solution that, normally, remains primal and dual feasible if the basis is inverted and used to solve for the primal and dual solution.

If crossover isn't run, then the primal and dual solution - because it comes from interior point - is far from guaranteed to satisfy primal and dual tolerances, and have small residual errors in the primal and dual equations. I found this out in the autumn when I added the code that determines a primal-dual solution that satisfies complementarity. I just accept some nontrivial residual errors. The lack of (exact) complementarity and residual errors in the primal-dual solution from interior point is what makes crossover such a challenge in general. IPX gets a more accurate solution from interior point than a Cholesky-based solver due to Lukas' smart preconditioning, so IPX crossover is usually trouble-tree.

[lgottwald]

That makes all sense, I understand now why that is how it is. For me a workaround on the MIP side is fine for the moment, but I would put it somewhere else than in the normal run() function. I'll let you know when I finished this.

[lgottwald]

For the basis from IPX I would, however, also add that it is set to not be alien if that is ok?