How do Quadratic Programming solvers handle variables without bounds?

Solvers for non-convex QPs generally do the McCormick relaxation of the term $$xy=z$$ and then do spatial branch and bound.

This requires that $$x$$ and $$y$$ have given bounds, how do they handle the case when the variables are unbounded?

Unless we can derive bounds during presolving, the standard way is to set a default variable range instead (e.g. $$\pm1.e16$$) so that we can generate the McCormick constraints.

There are numerical & convergence issues to consider depending on the number that is chosen, so setting smaller bounds might be preferable if we are certain that it's safe from a modelling point of view.

In general-purpose solvers we don't have that luxury because this is a modelling issue, so we pick large bounds by default and compensate for potentially bad numerics using a large bag of tricks, e.g., domain reduction techniques, re-scaling the problem, or higher precision arithmetic.

• 1e16? Ouch (see 1st sentence of @Johan Löfberg 's answer) , i hope that's in at least quad precision ... unless there is some special logic and handling to not get into trouble with that. Feb 18 '21 at 15:43
• @MarkL.Stone We use multiprecision arithmetic, but this actually works fine even with doubles 99% of the time because domain reduction techniques typically narrow the domain very quickly. There are many tricks to handle edge cases for bad numerics (too many to list here), but these are things we do in general for box-constrained problems, not just for bilinear relaxations. Feb 18 '21 at 17:24
• Thanks for adding the important clalrifying remarks. Feb 18 '21 at 17:56

Besides simply adding a large bound (which can cause numerical issues and lead to poor branching) or presolve from constraints involving the unbounded variable, the solver might be able to derive bounds by bound-propagation once a feasible solution is available.

As an example (assuming you have some bound, otherwise an indefinite problem will be unbounded), if you want to minimize $$-xy + y^2$$ and only have the bound $$0\leq x \leq 1$$ and you find the solution $$x = y = 0$$, you immediately have that the optimal solution satisfies $$-xy+y^2\leq 0$$ or equivalently $$y^2\leq xy$$. If $$x=0$$ it must hold that $$y=0$$, otherwise $$y$$ is non-negative, and by maximizing over $$x$$ to upper bound the right-hand side it holds that $$y^2 \leq y$$ which means $$y\leq 1$$, so in summary $$0 \leq y \leq 1$$ in any optimal solution.

I asked this question to the Gurobi support and this is their response:

Gurobi tries to deduce valid bounds for variables in presolve. Depending on the bounds of the participating variables, it might be possible to construct a partial McCormick relaxation.

If Gurobi can't find any bounds for one of the variables, the unboundedness is handled via branching. For example, Gurobi could branch on the unbounded variable at a large value like 1e6. The >=1e6 node is deferred for later - this is a case where you would see a "postponed" message in the log file. If Gurobi finds a feasible point, it can try to derive bounds for the deferred node through constraint propagation. It may also be able to compute a certificate of infeasibility for the node. Otherwise, the solver doesn't have much choice but to keep branching.