ERGO-Code
Strange bug in solving a QP model
Hi everyone.
First of all, thanks for your amazing work with this Solver.
I'm currently using HiGHS in my activities and recently I faced a strange beavhiour with a specific model. The model is attached to this message and it is in a txt format only because GitHub doen't allows to attach .lp files. HiGHS_model_strange.txt
It is a QP model, with just one constraint and all the variables are bounded beetween [-2,2]. For this reasons we should expect a finite objective value and in fact other solvers, like CPLEX, produce an optimal solution. However, you can verify that if you read and then run the model, HiGHS actually return an Unbounded model status. Can someone help me with this problem?
Moreover, if we slightly change something in the model, like some coefficients in the constraint, the model is actually solved to an optimal solution. I also attach this version of the model to prove this beavhiour. HiGHS_model_normal.txt
Thanks in advance for your time!
Enrico Calandrini
Running this (model strange) on my machine:
Running HiGHS 1.5.3 [date: 2023-07-27, git hash: ccb115d11-dirty] Copyright (c) 2023 HiGHS under MIT licence terms QP issue1363 has 1 rows; 18 cols; 9 matrix nonzeros; 18 Hessian nonzeros Iteration, Runtime, ObjVal, NullspaceDim 0, 0.000782, 4.542738, 0 28, 0.003547, -6.048694, 3 Model status : Optimal Objective value : 6.0486942966e+00 HiGHS run time : 0.00
How are you running the model? via API, or directly via command line?
Sorry, I forgot to specify. I am using the C API, thus the model is simply loaded and run by using the file attached: try_model.txt
What I get running this simple file is:
Running HiGHS 1.5.3 [date: 2023-07-11, git hash: 9e519472f] Copyright (c) 2023 HiGHS under MIT licence terms Iteration, Runtime, ObjVal, NullspaceDim 0, 0.001211, 4.542738, 0 Model status : Unbounded QP ASM iterations: 10 Objective value : -nan HiGHS run time : 0.00 Model status is 10
Thanks again.
This is an interesting one. At iteration 6, a (minor-iteration) search direction is computed that's not a descent direction. At iteration 7, a -NaN value appears in the search direction leading to the unboundedness-observation. I suspect the two issues are connected.
I have yet to find out whether that's an algorithmic problem (potentially in the pricing code, or the code handling semi-definiteness), or a numerical issue.
switching to a different pricing strategy helped in solving that particular instance, but that's likely down to luck.
Thank you very much for your work!
Let me know if I could help in any way. I'll keep an eye on this issue to see if you'll manage to resolve it.
I've spent a lot of time on this. Insights so far: The first couple of iterations are essentially flipping bounds for some active bound constraints. The basis consists of the constraint as well as all bounds except for the 9th variable. At iteration 5, bound 17 is dropped from active set without a new constraint becoming active. The reduced hessian is positive definite for this situation. at iteration 6, bound 8 is dropped from the active set and bound 9 enters the basis. either 17 or 8 can be chosen to leave the basis (both are in the "inactive" part of the basis. no matter which one you chose, the new nullspace after updating the basis leads to the reduced hessian being [0]. I am not sure if that is a numerical issues, a bug, or an algorithmic issue.
the reduced hessian being 0 then results in a -inf and a NaN value a bit later, which in turn leads to the false unboundedness statement by the solver.
I redid the basis calculations by hand (not trivial with an 18x18 matrix), and I just can't see what's wrong.
I am strongly suspecting that this error is the reason for a number of other positive-semidefinite QP instances failing, and this report is the smallest instance to debug. However, for the moment I spent way more time on this than I should have. I'll need to do other work for a while, but I will revisit this in the future.
I did find a number of other issues while debugging this, and I'll create a pull request to fix those.
If you don't mind explaining, I'd be interested in learning about the context of this instance. It looks like a computer-generated subproblem.
Thank you very much for your work!
Of course! At the moment I'm collaborating in the development of a modeling system named SMS++. In particular, I have worked in building the interface to allow users to solve models using HiGHS (here you can see that we developed various *MILPSolver using different commercial and open-source solvers). As a matter of fact, if you have time and are interested in this, you could give a look to HiGHSMILPSolver.cpp and give me your feedback! I would really appreciate this.
More in detail, this instance has been generated by our tester LagrangianDualSolver_Box. Essentially it creates a random model (LP or QP) with a number of constraints specified by the user. This problem is then repeatedly randomly modified and re-solved several times. You can find a detailed description of this tester in the link above.
Of course! At the moment I'm collaborating in the development of a modeling system named SMS++. In particular, I have worked in building the interface to allow users to solve models using HiGHS (here you can see that we developed various *MILPSolver using different commercial and open-source solvers). As a matter of fact, if you have time and are interested in this, you could give a look to HiGHSMILPSolver.cpp and give me your feedback! I would really appreciate this.
Ah, I'm familiar with Antonio and SMS++. I've had a quick look at HiGHSMILPSolver.cpp, and it looks like a well written interface. However, curious as to whether you have the functionality to pass a (partial) integer feasible solution to HiGHS, I searched for setSolution, and spotted that the only occurrences are the following.
if( Highs_setSolution( highs , NULL , NULL , dj.data() , NULL ) == kHighsStatusError ) throw( std::runtime_error( "Unable to get reduced costs with Highs_setSolution") );
I note that the call should be to Highs_getSolution
Hence I deduce that you don't have the functionality to pass a (partial) integer feasible solution to HiGHS. If you do have a heuristic to generate such a solution, it could reduce MIP solution time significantly.
More in detail, this instance has been generated by our tester LagrangianDualSolver_Box. Essentially it creates a random model (LP or QP) with a number of constraints specified by the user. This problem is then repeatedly randomly modified and re-solved several times. You can find a detailed description of this tester in the link above.
Experience makes me question the validity of performing experiments with random problems in anything other than simple situations of continuous optimization. If you're only testing the logical correctness of code, then it's OK. Otherwise, in what sense are you testing the system?
Hi, sorry for the delay but I have talked to prof. Frangioni about your suggestions and we are seriously thinking about adding the possibility to pass integer feasible solution to HiGHS. It could be a nice improvement!
I've had a quick look at HiGHSMILPSolver.cpp, and it looks like a well written interface.
Thank you very much!
If you're only testing the logical correctness of code, then it's OK.
Yes thats exactly what we are using those tests for. In SMS++ we are focusing the attention on solving not a single instance of a problem but a lot of different modification of a single problem, verifying that in each step solvers are actually able to give the correct answer.