I want to solve a mixed integer quadratic programming problem with 267 variables [1] using SCIP.
CPLEX can solve the problem in about 30 seconds and a solution that is extremely close to the optimum is already found within a fraction of a second [2, 3].
Unfortunately, SCIP really struggles with this problem, unable to find a solution that is anywhere near the optimum even after running for over 20 minutes [4].
Why is this? Is CPLEX really that much better at MIQP than SCIP? Did I not configure SCIP correctly? How can I solve this problem with SCIP?
It also seems to me like the solutions SCIP finds are very far away from the solution of the relaxation. I was under the impression that SCIP would first solve the relaxation and then try to find an integer solution based on that. Is this not correct? If yes, why are the solutions so far off?
Why is this?
Your problem has a quite dense convex quadratic objective function in all variables, one linear constraint, and all variables are restricted to be integral.
SCIP recognizes the convexity, but then tries to construct a linear relaxation by alternating between solving the LP and adding a linearization of the quadratic function to the LP. It uses the solution of the LP as reference point (like in Kelley's cutting plane method). This can mean slow convergence (and better approaches are to linearize at the boundary of the feasible region), that is, requiring many LP solves until you get to a decent bound. In addition, the linearization cuts that are generated are dense, which means that with every new cut, the LP gets slower and slower. After a number of rounds in the root node, SCIP decides that separation does not make enough progress (stalling) and it goes into branching on fractional integer variables. However, since the relaxation is still so bad, the dual bound doesn't give any good guidance on which variable to branch on. And since the LP solve is so slow, it doesn't get to process nodes quickly.
CPLEX very likely solves the convex QP relaxation of the MIQP to compute a dual bound. This gives a very good bound close to the optimal value. I don't know if it then uses this point to linearize the quadratic function, which would result in the LP bound to be equal to the QP bound, and then continues with the LP, or it continues by solving the QP.
Is CPLEX really that much better at MIQP than SCIP?
On this instance (and probably other convex MIQP with dense objective), yes.
Did I not configure SCIP correctly?
No.
How can I solve this problem with SCIP?
There is no magic switch in SCIP to just make it behave like CPLEX.
Switching from the LP relaxation to the QP relaxation is something that is not available in the main SCIP. Something that comes close is the NLP diving heuristic. This solves the QP relaxation, then fixes or bounds some integer variable with fractional value, and then solves again. It needs quite a deep dive at the root node to get to a feasible solution, but since there are not many constraints, it is possible and I get to a solution with value 52474.45. Unfortunately, in default settings the heuristic does not run early enough. If you set
heuristics/nlpdiving/freq = 0
heuristics/nlpdiving/freqofs = 0
heuristics/nlpdiving/maxnlpiterabs = 10000
heuristics/nlpdiving/nlpfastfail = FALSE
heuristics/nlpdiving/varselrule = f
it will run once at the end of the root node (to make it run earlier, change HEUR_TIMING in heur_nlpdiving.c to SCIP_HEURTIMING_BEFORENODE).
This still leaves the problem that the LP relaxation is useless. So SCIP would not be able to proof optimality in a reasonable amount of time. You can get a better bound if you also set
constraints/nonlinear/linearizeheursol = e
With this, a linearization of the quadratic function in the solution found by the NLP diving heuristic is added to the LP. With that, I get to a gap of 8% at the end of the root node.
Log: https://gist.github.com/svigerske/dfdb9e95af1f4eb386b8b05770ffde4c
While there is no QP relaxation in the main SCIP, there is an example that shows how to use the NLP relaxation for bounding. You can find this under examples/Relaxator. With this, you can get closer to what CPLEX is doing, but this is indeed just a very small example.
If one switches off the LP relaxation and uses the NLP relaxation, you get a good dual bound in the root node (52467.99), but no primal bound anymore. This is because many primal heuristic work with the LP relaxation and are essentially off now. As before one can enable the NLP diving heuristic to get a good primal bound (52597.73). However, closing the remaining 0.25% gap seems difficult. First, solving the NLP relaxation is not fast either, because we do not (and cannot really) warmstart Ipopt. Second, there is no sophisticated branching rule implemented.
If you want to try out the NLP relaxator example for your instance, do these changes to the source:
--- a/examples/Relaxator/src/relax_nlp.c
+++ b/examples/Relaxator/src/relax_nlp.c
@@ -32,8 +32,8 @@
#define RELAX_FREQ 1
#define NLPITERLIMIT 500 /**< iteration limit of NLP solver */
-#define FEASTOLFAC 0.01 /**< factor for NLP feasibility tolerance */
-#define RELOBJTOLFAC 0.01 /**< factor for NLP relative objective tolerance */
+#define FEASTOLFAC 1.0 /**< factor for NLP feasibility tolerance */
+#define RELOBJTOLFAC 1.0 /**< factor for NLP relative objective tolerance */
/*
* Data structures
diff --git a/src/scip/heur_nlpdiving.c b/src/scip/heur_nlpdiving.c
index 88a9237d12..bc1ab8f690 100644
--- a/src/scip/heur_nlpdiving.c
+++ b/src/scip/heur_nlpdiving.c
@@ -63,7 +63,7 @@
#define HEUR_FREQ 10
#define HEUR_FREQOFS 3
#define HEUR_MAXDEPTH -1
-#define HEUR_TIMING SCIP_HEURTIMING_AFTERLPPLUNGE
+#define HEUR_TIMING SCIP_HEURTIMING_BEFORENODE
#define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
/* event handler properties */
and set the options
lp/solvefreq = -1
relaxing/nlp/priority = 0
relaxing/lp/freq = -1
heuristics/nlpdiving/freq = 0
heuristics/nlpdiving/freqofs = 0
heuristics/nlpdiving/maxnlpiterabs = 10000
heuristics/nlpdiving/nlpfastfail = FALSE
Log: https://gist.github.com/svigerske/5ebb29c5fc6425a60b023944f5b07c04
I was under the impression that SCIP would first solve the relaxation and then try to find an integer solution based on that. Is this not correct? If yes, why are the solutions so far off?
If the relaxation is far off, then it doesn't give a good guidance on finding an integer solution nearby. The LP relaxation consists of your linear constraints, the bounds on the integer variables, and some cuts that very badly approximate the quadratic objective function. The primal heuristic that use the LP relaxation find a feasible solution for this, but the objective function value of the relaxation is at around -4e6. There is then something that updates the objective function value to use the original quadratic function, but this ends up at around 2e6, far from any optimal value.
Other free/open-source solvers that may handle this instance well and come into my mind are SHOT and MINOTAUR. The former should compute tight linear relaxations, but I don't know how well it does on the primal bound. The latter is specialized on solving QP relaxations quickly and do QP diving.