1. The problem
Draft
Four letters
This subsection fixes the objects every later section is about: the problem, its feasible set and optimal value, the four classes the letters name, and the one device, relaxation, on which every bound in the series rests. The definitions are elementary. They are stated with care because the words "optimal", "bound" and "gap" carry precise meanings in a solver's output, and a reader who takes them loosely will misread that output.
Definition 1.1.1 (Optimization problem, feasible set, optimum). An optimization problem consists of variables, constraints that the variables must satisfy, and an objective function to be made as small as possible. In this series the variables are a vector \((x, y) \in \mathbb{R}^n \times \mathbb{Z}^p\), the constraints are inequalities \(g_i(x, y) \le 0\), \(i = 1, \dots, m\), together with bounds \(l \le (x, y) \le u\), and the feasible set \(\mathcal F\) is the set of all \((x, y)\) that satisfy every constraint. The optimal value is
\[z^\star \;=\; \inf\,\{\, f(x, y) : (x, y) \in \mathcal F \,\}.\]An optimal solution, or global minimizer, is a feasible point at which the infimum is attained. The problem is infeasible if \(\mathcal F\) is empty, in which case \(z^\star = +\infty\), and unbounded if \(z^\star = -\infty\). An equality constraint \(h(x, y) = 0\) is the pair of inequalities \(h \le 0\) and \(-h \le 0\), and maximizing \(f\) is minimizing \(-f\), so the display loses nothing. When the feasible set is nonempty, every variable has finite bounds and the functions are continuous, the feasible set is a closed subset of a box, hence compact, and the infimum is attained by the Weierstrass theorem. Without finite bounds on the integer variables the infimum of even a convex problem can fail to be attained, and Section 1.5 gives the standard example. This is the first of several reasons that global solvers insist on bounds.
Definition 1.1.2 (The four classes). The problem of Definition 1.1.1 is a linear program (LP) if \(p = 0\) and \(f\) and every \(g_i\) are affine. It is a mixed-integer linear program (MILP) if \(f\) and every \(g_i\) are affine and \(p \ge 1\). It is a nonlinear program (NLP) if \(p = 0\) and some function is not affine, and a mixed-integer nonlinear program (MINLP) if \(p \ge 1\) and some function is not affine. An NLP or MINLP is convex if \(f\) and every \(g_i\) are convex functions on the box \(B = [l, u]\) in the sense of Definition 1.2.1, so that dropping the integrality of \(y\) leaves a convex program. Otherwise it is nonconvex. Two conventions are worth stating. An integer variable with bounds \(0\) and \(1\) is a binary variable, and most of the integer variables in this series are binary, because they encode decisions. A MINLP whose nonlinear functions are all quadratic is a mixed-integer quadratically constrained program (MIQCP), and a MINLP whose only nonlinearity is a quadratic objective is a mixed-integer quadratic program (MIQP). Both appear in Sections 4 and 9.
The table sorts the four classes by what is known about solving them. It uses one word that needs fixing first. A certificate is a finite piece of data, such as a vector of multipliers or a dual solution, from which a claim about the problem can be checked without solving the problem again. The table's two right-hand cells are the subject of the series.
| linear functions | nonlinear functions | |
|---|---|---|
| continuous variables | LP: exact; polynomial in theory (ellipsoid, barrier), and the simplex method, exponential in the worst case, fast in practice | NLP: local methods find a point satisfying the first-order (KKT) conditions, a candidate for a local minimizer; local optimality is certified only under second-order conditions, global optimality needs convexity or a search |
| some integer variables | MILP: branch and bound on LP relaxations, plus cuts; NP-hard, routinely solved | MINLP: branch and bound on NLP or on a convex relaxation; both difficulties at once; routinely not solved |
The LP cell says three things, and the order matters. Linear programming is solvable in polynomial time, by the ellipsoid method and by interior-point (barrier) methods. The simplex method, which is what most solvers run on the relaxations of a tree, is exponential in the worst case and fast in practice. Smoothed analysis explains the second half: the expected running time of the simplex method with the shadow-vertex pivot rule is polynomial under small random perturbations of any input.L. G. Khachiyan, "Polynomial algorithms in linear programming", USSR Computational Mathematics and Mathematical Physics 20 (1980), the full account of the ellipsoid method, announced in a 1979 note in Doklady Akademii Nauk SSSR; N. Karmarkar, "A new polynomial-time algorithm for linear programming", Combinatorica 4 (1984); D. A. Spielman and S.-H. Teng, "Smoothed analysis of algorithms: why the simplex algorithm usually takes polynomial time", Journal of the ACM 51 (2004). An LP's answer is exact, in the sense that for rational data an optimal vertex has rational coordinates of polynomial size and a solver can, in principle, return it exactly. What a floating-point solver returns in practice is treated in Sections 5.6 and 7.3.
The NLP cell is the one most often misread. A local method for a nonlinear program, such as an interior-point or sequential quadratic programming code, stops at a point that satisfies the first-order optimality conditions of Karush, Kuhn and Tucker to a tolerance. That point is a candidate for a local minimizer and nothing more. In a nonconvex problem it may be a saddle point, or a local minimizer far above the global one.S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004), chapter 5, for the conditions and for the fact that they are sufficient under convexity. Section 1.3 states the conditions, gives a two-variable example with a KKT point that is a saddle, and says what can be certified about a local solver's answer. When the problem is convex the same point is a global minimizer, and this is the whole reason the convex case is tractable. The MILP cell contains the first of the two difficulties the series is about. The problem is NP-hard, so no known method avoids, on its worst instances, a number of steps exponential in the number of integer variables. The reason is that linear inequalities over \(0\)–\(1\) variables can say anything a Boolean formula can say. The clause "\(x_1\) or not \(x_2\) or \(x_3\)" is the inequality \(x_1 + (1 - x_2) + x_3 \ge 1\), one such inequality per clause turns any formula into a system, and the formula is satisfiable exactly when the system has a \(0\)–\(1\) point. A method that decided the feasibility of every such system in polynomial time would therefore decide satisfiability in polynomial time, which is the problem that no one has done and that P \(\ne\) NP forbids. Theorem 1.5.2 states this precisely, and Section 1.5 says what the statement does and does not imply: it is about the worst case over all instances of growing size, and it says nothing about any particular instance. Instances with tens of thousands of such variables are nevertheless solved every day, by branch and bound on linear relaxations strengthened with cutting planes, linear inequalities that every integer point satisfies but the relaxation's current solution violates (Section 3.3).R. M. Karp, "Reducibility among combinatorial problems", in Complexity of Computer Computations (Plenum, 1972), where 0–1 integer programming is one of the twenty-one problems; Section 1.5 says precisely what NP-hardness does and does not imply. Branch and bound is A. H. Land and A. G. Doig, "An automatic method of solving discrete programming problems", Econometrica 28 (1960), and the two-child form every solver uses is R. J. Dakin, "A tree-search algorithm for mixed integer programming problems", The Computer Journal 8 (1965); Section 3.1 gives the algorithm. How far the linear solvers have come is measured in Koch, Berthold, Pedersen and Vanaret (2022), cited in the introduction and reported in Section 5.1. The MINLP cell inherits that difficulty and adds the second.
(The certificate, cell by cell) The second difficulty is best seen by asking what each cell's certificate actually is. Branch and bound, the search of Section 3.1, splits the problem into subproblems, its nodes, computes a bound for each, and prunes a node, that is, drops it from the search, when the node's bound is no better than the incumbent, the best feasible point found so far. So the number computed at a node has to be a bound, a value the optimum over that node cannot beat (Proposition 1.1.4 below), and the solver has to know that it is one. In the LP cell the knowledge is a dual solution: for the program \(\min\{c^\top x : Ax \ge b,\ x \ge 0\}\), a vector \(y \ge 0\) with \(A^\top y \le c\) proves that \(b^\top y \le c^\top x\) for every feasible \(x\) (Theorem 2.2.6), and checking this is a matrix-vector product, not a second solve. The relaxation at every node of a MILP tree is an LP, so every node's bound comes with such a vector, which is why the MILP cell has only the first difficulty. In the convex half of the NLP cell the certificate is the vector of KKT multipliers at the point the local method returns, because under convexity the KKT conditions are sufficient (Proposition 1.3.5), so a point and multipliers that satisfy them prove the point is a global minimizer, and dropping the integrality of a convex MINLP leaves a convex program whose local solution is therefore a certified bound. Branch and bound then runs on NLP relaxations exactly as it runs on LP relaxations, and this is the convex MINLP of Section 3.4: hard for the reason MILP is hard, and for no other.
(Why a local method's answer is not a bound) When the nonlinear functions are not convex, the relaxation obtained by dropping integrality is a nonconvex NLP, and nothing a local method returns from it is a certificate of a bound. The method stops at a KKT point, and the multipliers there prove only that the point is stationary. Its objective value is the value of one feasible point of the relaxation, which is an upper bound on the relaxation's optimum and a lower bound on nothing; and the method may stop at an infeasible point instead, since in a nonconvex feasible set a local method cannot certify infeasibility either. A tree that pruned on such a value would discard subtrees containing the optimum, with no sign in its output that it had done so. The point still has a use. If it satisfies the integrality constraints it is an incumbent, and that is how every solver of Section 5 uses local methods, as suppliers of candidates for the upper side of the bracket. On the lower side it is of no use at all.
(Relaxing twice, and branching on a continuous variable) A global solver therefore relaxes twice. It drops integrality, and it replaces every nonconvex function by a convex function that lies beneath it on the current box, so that the problem solved at the node is convex and its multipliers certify its value as a lower bound. Section 2.4 constructs these underestimators term by term; for the product \(xy\) on a rectangle they are the four McCormick planes that Section 1.6 applies to the running example R2. The price is that an underestimator is a function of the box. On a wide box it can lie far beneath the function it replaces, and the only way to bring it closer is to make the box smaller, which no branching on integer variables can do. So the search splits the range of a continuous variable, relaxes again on each half, and continues until the gap between the bound and the incumbent is below tolerance. This is the spatial branch and bound of Section 3.5, and it is why the table's MINLP cell says "on NLP or on a convex relaxation": on the NLP when the functions are convex, on a convex relaxation of the NLP when they are not.Branch and bound on continuous variables begins with J. E. Falk and R. M. Soland, "An algorithm for separable nonconvex programming problems", Management Science 15 (1969), for separable functions, and takes its general form in McCormick's 1976 paper, cited in Section 1.6, which supplies the underestimators. R. Horst and H. Tuy, Global Optimization: Deterministic Approaches (Springer, 3rd ed., 1996) is the reference for the convergence theory, where the condition on the splitting is that it be exhaustive, that is, that the boxes along any infinite path of the tree shrink to a point.
(Finite only up to a tolerance) Branching on a continuous variable changes what finiteness means. An integer variable with finite bounds can be branched on only finitely often, since each branch removes at least one integer from its range, so a MILP tree is finite and the algorithm ends with the exact optimum. A continuous range can be halved without end. The tree a global solver builds is finite only because the solver stops splitting a node once its bound is within \(\varepsilon\) of the incumbent, or once the box is narrower than a tolerance, and what it returns is a point and a bound at most \(\varepsilon\) apart, not the optimum itself. Section 1.5 states this certificate exactly, as Definition 1.5.20, and it is the form every result reported in Section 5 takes. A nonconvex NLP with no integer variables already needs all of this. The pooling problem R4 of Section 1.6, whose only nonconvexity is two products of continuous variables, has nothing to branch on but a continuous variable and belongs in the MINLP cell for every practical purpose; a true MINLP adds the first difficulty on top of it.
The device that connects the four cells is relaxation. The series defines it fully in Section 2.1 and uses it everywhere, so a preview belongs here.
Definition 1.1.3 (Relaxation). A relaxation of the minimization problem \(\min\{f(x) : x \in \mathcal F\}\) is a second minimization problem \(\min\{f_R(x) : x \in \mathcal F_R\}\) with \(\mathcal F \subseteq \mathcal F_R\) and \(f_R(x) \le f(x)\) for every \(x \in \mathcal F\). Its optimal value is written \(z_R\).
Proposition 1.1.4 (A relaxation bounds the optimum). If the second problem is a relaxation of the first, then \(z_R \le z^\star\). If moreover a minimizer \(\bar x\) of the relaxation lies in \(\mathcal F\) and \(f_R(\bar x) = f(\bar x)\), then \(\bar x\) is a global minimizer of the original problem and \(z_R = z^\star\).
Proof. Every \(x \in \mathcal F\) lies in \(\mathcal F_R\), so \(z_R \le f_R(x) \le f(x)\). Taking the infimum over \(x \in \mathcal F\) gives \(z_R \le z^\star\). For the second statement, \(f(\bar x) = f_R(\bar x) = z_R \le z^\star \le f(\bar x)\), so every inequality is an equality. ∎
The proposition is all that a bound from a relaxation asserts: a number computed from an easier problem that the true optimum cannot beat. Section 2.2 explains why such a number is called a dual bound. The LP relaxation of a MILP, obtained by deleting \(y \in \mathbb{Z}^p\), is a relaxation with \(f_R = f\). The convex relaxations of Section 2 replace a nonconvex \(g_i\) or \(f\) by a convex function lying beneath it, which enlarges \(\mathcal F\) and lowers \(f\) at once. The second half of the proposition is the termination test of every exact method. When the relaxation's answer happens to be feasible for the original problem, the search at that node is over.
Since this is the first time the definitions are applied, it is worth writing out what each symbol is on the running example R3 of Section 1.6, the MINLPLib instance st_e13,
\[\min\ 2x + y \quad \text{subject to}\quad 1.25 - x^2 - y \le 0,\quad x + y \le 1.6,\quad 0 \le x \le 1.6,\quad y \in \{0, 1\}.\]In the notation of Definition 1.1.1, \(n = p = 1\), the variables are \((x, y) \in \mathbb{R} \times \mathbb{Z}\), the objective is \(f(x, y) = 2x + y\), the constraints are \(g_1(x, y) = 1.25 - x^2 - y\) and \(g_2(x, y) = x + y - 1.6\), and the bounds are \(l = (0, 0)\) and \(u = (1.6, 1)\). Since \(p \ge 1\) and \(g_1\) is not affine, the problem is a MINLP by Definition 1.1.2, and since \(g_1\) is concave in \(x\) it is a nonconvex one. The feasible set is what the three conditions leave: \(y = 1\) forces \(x \ge 0.5\) from \(g_1\) and \(x \le 0.6\) from \(g_2\), while \(y = 0\) forces \(x \ge \sqrt{1.25}\), so
\[\mathcal F \;=\; \{(x, 1) : 0.5 \le x \le 0.6\} \;\cup\; \{(x, 0) : \sqrt{1.25} \le x \le 1.6\},\]two segments. The objective is smallest at the left end of each, \(2.0\) at \((0.5, 1)\) and \(2\sqrt{1.25} = 2.236\) at \((1.118, 0)\), so \(z^\star = 2.0\), attained at the global minimizer \((0.5, 1)\). The simplest relaxation of R3, the one Section 2.4 will construct in general, is the chord relaxation. On a box \(0 \le x \le u\) the term \(-x^2\) is replaced by its chord, the line \(-ux\) through \((0, 0)\) and \((u, -u^2)\), which lies beneath the parabola on the box because \(x^2 \le ux\) there; and the integrality of \(y\) is dropped. In the notation of Definition 1.1.3,
\[f_R = f, \qquad \mathcal F_R \;=\; \{(x, y) : 1.25 - ux - y \le 0,\ x + y \le 1.6,\ 0 \le x \le u,\ 0 \le y \le 1\},\]and \(\mathcal F \cap \{x \le u\} \subseteq \mathcal F_R\) because \(1.25 - ux - y \le 1.25 - x^2 - y\) on the box, so every point that satisfies \(g_1 \le 0\) satisfies the relaxed constraint too. The relaxed problem is a linear program, and its minimum is on the line \(y = 1.25 - ux\) at the largest \(y\) allowed, \(y = 1\), giving \(z_R = 1 + 0.5/u\) at \((0.25/u, 1)\): at the root box \(u = 1.6\) this is \(1.3125 \le 2.0 = z^\star\), the first half of the proposition, and at \(u = 0.5\) the relaxed minimizer is \((0.5, 1) \in \mathcal F\) with \(f_R = f\), the second half.
Relaxation connects the cells: deleting y in Z^p moves a problem up
linear functions nonlinear functions
continuous LP NLP
variables ^ ^
| |
| delete y in Z^p: | delete y in Z^p: a
| the LP relaxation, | convex program when f
| f_R = f | and every g_i are
| | convex; otherwise no
| | local method solves it
| | to a certificate
some integer MILP MINLP
variables
One convention concerns signs. Every general display in this series minimizes, so a relaxation's value is a lower bound and an incumbent's value is an upper bound. Three running examples, the drawn two-variable program R1, the bilinear example R2 and the pooling problem R4, are written as maximizations in Section 1.6, because that is how they are naturally stated, and for them the inequalities reverse. The text says so each time one of them appears, and Section 2.1 states the reversed chain of bounds for a maximization, so that the figures' readouts can be read against the text.
The proposition gives one side of a bracket, and the other side needs one more word. An incumbent is the best feasible point a solver has found so far, and its objective value, written \(z_{\mathrm{inc}}\), is an upper bound on \(z^\star\) for no deeper reason than that \(z^\star\) is the infimum of \(f\) over \(\mathcal F\) and the incumbent is one point of \(\mathcal F\). So at every moment of a search the optimum is pinned between two numbers, \(z_R \le z^\star \le z_{\mathrm{inc}}\), the first computed from a relaxation and the second read off a feasible point, and the difference \(z_{\mathrm{inc}} - z_R\) is the gap. A solver works on both ends at once: it raises \(z_R\) by branching and by tighter relaxations, and it lowers \(z_{\mathrm{inc}}\) by finding better feasible points, and it stops when the two meet, or come within a tolerance of meeting. Section 2.1 defines the gap and its relative version precisely, and Section 2.3 reads a solver's log in these terms.
Proposition 1.1.4 on the objective axis, in both senses
minimizing (every general display): z_R <= z* <= z_inc
------+--------------+------------------------+--------->
z_R z* z_inc
the relaxation: the incumbent:
a lower bound an upper bound
maximizing (R1, R2 and R4, Section 1.6): the chain reverses,
z_inc <= z* <= z_R
------+------------------------+--------------+--------->
z_inc z* z_R
the incumbent: the relaxation:
a lower bound an upper bound
a minimizer xbar of the relaxation that lies in F, with
f_R(xbar) = f(xbar): z_R = z* = f(xbar), and the search at that
node is over
Where this is used
A solver's log prints the best bound and the incumbent at every step. The other columns describe how they were obtained. Section 2.1 defines the gap between them and the conventions by which different solvers report it, and Section 2.3 walks through one complete solve, line by line, on the drawn example R1.
Convexity, not linearity
The class table of the previous subsection is organized by linearity, and that is not where the difficulty lies. This subsection draws the line that does organize the subject, between convex and nonconvex problems. It states the theorem that makes the convex side easy, and it shows on two sets described by one equation how a curve can fall on either side of the line. The figure at the end puts the two halves of the statement side by side.
Definition 1.2.1 (Convex set, convex function, epigraph). A set \(S \subseteq \mathbb{R}^n\) is convex if for every \(x, y \in S\) and every \(t \in [0, 1]\) the point \((1 - t)x + ty\) lies in \(S\): the segment between any two points of \(S\) stays in \(S\). A function \(f : S \to \mathbb{R}\) on a convex set is convex if
\[f\big((1 - t)x + ty\big) \;\le\; (1 - t)\, f(x) + t\, f(y) \qquad \text{for all } x, y \in S,\ t \in [0, 1],\]strictly convex if the inequality is strict for \(x \ne y\) and \(t \in (0, 1)\), and concave if \(-f\) is convex. The epigraph of \(f\) is \(\operatorname{epi} f = \{(x, s) : x \in S,\ s \ge f(x)\}\), and \(f\) is convex if and only if \(\operatorname{epi} f\) is a convex set. The convex hull \(\operatorname{conv}(T)\) of an arbitrary set \(T\) is the smallest convex set containing it, equivalently the set of all finite convex combinations \(\sum_k t_k x_k\) with \(x_k \in T\), \(t_k \ge 0\), \(\sum_k t_k = 1\).
The equivalence between a convex function and a convex epigraph is the reason the two notions can be treated as one. The inequality says that the chord of the graph of \(f\) over \([x, y]\) lies above the graph. The graph lies below all of its chords exactly when the epigraph contains the segment between any two of its points, which is the convexity of the epigraph. For a twice differentiable function on an open convex set, convexity is the same as a positive semidefinite Hessian, the matrix of second partial derivatives, everywhere. That is the test solvers apply to quadratic functions, and it is the reason the quadratic case has such a sharp theory.
Definition 1.2.2 (Convex problem, local and global minimizers). The problem \(\min\{f(x) : x \in \mathcal F\}\) is convex if \(\mathcal F\) is a convex set and \(f\) is a convex function on it. A point \(\bar x \in \mathcal F\) is a global minimizer if \(f(\bar x) \le f(x)\) for all \(x \in \mathcal F\), and a local minimizer if there is a neighbourhood \(U\) of \(\bar x\) with \(f(\bar x) \le f(x)\) for all \(x \in \mathcal F \cap U\). When \(\mathcal F\) is described by constraints \(g_i(x) \le 0\) with every \(g_i\) convex, it is a convex set (Proposition 1.2.5). That is how the convexity of a problem is recognized in practice: through its functions, not through its feasible set directly.
Theorem 1.2.3 (For a convex problem every local minimizer is global). Let \(\mathcal F\) be convex and \(f\) convex on \(\mathcal F\). Then every local minimizer of \(f\) over \(\mathcal F\) is a global minimizer, the set of global minimizers is convex, and if \(f\) is strictly convex the global minimizer is unique when it exists.The theorem and the two propositions that follow it are classical, and the proofs given here are the standard ones. Boyd and Vandenberghe (2004), cited above: section 4.2.2 for Theorem 1.2.3; section 3.1.3 for the first-order condition of Proposition 1.2.4; section 3.1.6 for the sublevel sets and section 2.2.3 for the second-order cone of Proposition 1.2.5; exercise 4.26 for the identity between a hyperbolic constraint and a second-order cone constraint.
Proof. Let \(\bar x\) be a local minimizer with neighbourhood \(U\), and suppose some \(x \in \mathcal F\) has \(f(x) < f(\bar x)\). For \(t \in (0, 1]\) the point \(x_t = \bar x + t(x - \bar x)\) lies in \(\mathcal F\) by convexity of the set, and by convexity of the function
\[f(x_t) \;\le\; (1 - t)\, f(\bar x) + t\, f(x) \;<\; f(\bar x).\]As \(t \to 0\) the point \(x_t\) enters \(U\), which contradicts local minimality. So no such \(x\) exists and \(\bar x\) is global. If \(x\) and \(y\) are both global minimizers with value \(z^\star\), then every point of the segment has \(f \le (1 - t) z^\star + t z^\star = z^\star\), hence \(f = z^\star\) there, so the set of minimizers is convex. If \(f\) is strictly convex and \(x \ne y\) were two minimizers, the midpoint would have \(f < z^\star\), which is impossible. ∎
(Why the convex side has certificates) In a picture, a convex function has no valley other than the lowest one. From any feasible point the direction toward a better feasible point is one along which the function is lower than \(f(\bar x)\) at every point of the segment, so a method that only knows how to go downhill cannot be trapped. The theorem is the reason the convex side of the table has certificates. A local method on a convex problem stops at a local minimizer, which is global, and the stationarity conditions it verifies at that point are a proof. Section 1.3 gives the conditions and the proof. On the nonconvex side the same method stops at the first valley it reaches, and the figure below shows how often that is the wrong one.
Two further consequences of convexity are used constantly in Sections 2 and 3, and it is worth having them in hand now.
Proposition 1.2.4 (A convex function lies above its tangents). Let \(f\) be convex and differentiable on an open convex set \(S\). Then for all \(x, y \in S\),
\[f(y) \;\ge\; f(x) + \nabla f(x)^\top (y - x).\]Consequently, if \(g\) is convex and differentiable and \(\ell_x(y) = g(x) + \nabla g(x)^\top (y - x)\) denotes its linearization at \(x\), then every point with \(g(y) \le 0\) satisfies the linear inequality \(\ell_x(y) \le 0\), for every choice of \(x\). If \(g\) is concave the reverse inequality holds, so \(\{y : \ell_x(y) \le 0\} \subseteq \{y : g(y) \le 0\}\): the linearized constraint is a restriction rather than a relaxation. It removes every feasible \(y\) with \(\ell_x(y) > 0\), and for a strictly concave \(g\) these include every boundary point other than \(x\).
Proof. For \(t \in (0, 1]\), convexity gives \(f(x + t(y - x)) \le f(x) + t\,(f(y) - f(x))\), so \([f(x + t(y - x)) - f(x)]/t \le f(y) - f(x)\). Letting \(t \downarrow 0\), the left side tends to the directional derivative \(\nabla f(x)^\top (y - x)\). For the first consequence, \(\ell_x(y) \le g(y) \le 0\). For concave \(g\) apply the inequality to \(-g\): then \(g(y) \le \ell_x(y)\) for every \(y\), so \(\ell_x(y) \le 0\) implies \(g(y) \le 0\), which is the inclusion, and a feasible \(y\) is removed exactly when \(\ell_x(y) > 0\). If \(g\) is strictly concave then \(g(y) < \ell_x(y)\) for \(y \ne x\), so a boundary point \(y \ne x\), where \(g(y) = 0\), has \(\ell_x(y) > 0\). ∎
(Gradient cuts, and why they fail on concave constraints) The first consequence is a gradient cut. A cut, throughout the series, is a linear inequality that every feasible point satisfies, added to a relaxation to tighten it (Section 3.3), and a gradient cut is the tangent half-space of a convex constraint. A convex constraint implies each of its tangent half-spaces, so a convex feasible set is the intersection of the half-spaces of all its tangents and can be approximated from outside by finitely many of them. This is the engine of the outer-approximation methods of Section 3.4 and of every global solver that works with linear relaxations (Section 3.3). The second consequence is the reason those methods fail on nonconvex constraints. The tangent of a concave function lies above it, so the linearized inequality is a restriction rather than a relaxation, and a method that adds it as a cut may remove the optimum. Running example R3 in Section 1.6 is a two-variable instance on which exactly this happens, and Section 3.4 works it through. The arithmetic is short enough to show here. R3's curved constraint is \(g(x, y) = 1.25 - x^2 - y \le 0\), a concave \(g\), and its linearization at a point \((x_t, y_t)\) is \(\ell(x, y) = 1.25 + x_t^2 - 2x_t x - y\), so the inequality \(\ell \le 0\) reads \(y \ge 1.25 + x_t^2 - 2x_t x\). At the optimum \((0.5, 1)\) of R3 the right side is \(1.25 + x_t^2 - x_t = 1 + (x_t - 0.5)^2\), which exceeds \(1\) for every \(x_t \ne 0.5\). So every tangent but the one at the optimum itself cuts the optimum off, which is the second consequence with its boundary points made explicit.
Proposition 1.2.4 on a line: g against its linearization l_x at x
x a point with g(x) = 0; [ ] closed ends, < > unbounded
g convex g <= 0 [=================x
l_x <= 0 <=======================x
every feasible point is kept: a relaxation,
the gradient cut
g concave g <= 0 <=========] x==============>
l_x <= 0 x==============>
the feasible points with l_x > 0 are cut off: a
restriction; for g strictly concave they include
every boundary point other than x
Proposition 1.2.5 (Sublevel sets, and one hyperbola seen from two sides). (i) If \(g\) is convex on a convex set \(S\), then \(\{x \in S : g(x) \le 0\}\) is convex, and the intersection of convex sets is convex. Hence a problem whose constraint functions are all convex has a convex feasible set. (ii) In the open positive quadrant the set \(\{(x, y) : xy \ge 1\}\) is convex and the set \(\{(x, y) : xy \le 1\}\) is not. (iii) For \(x, y \ge 0\),
\[xy \;\ge\; 1 \quad\iff\quad \big\lVert (2,\ x - y) \big\rVert_2 \;\le\; x + y ,\]so the convex side of the hyperbola is a slice of a second-order cone.
Proof. (i) If \(g(x) \le 0\) and \(g(y) \le 0\) then \(g((1 - t)x + ty) \le (1 - t) g(x) + t\, g(y) \le 0\). An intersection of convex sets contains the segment between any two of its points because each set does. (ii) For \(x > 0\) the condition \(xy \ge 1\) reads \(y \ge 1/x\), and \(1/x\) is convex on \(x > 0\), since its second derivative is \(2/x^3 > 0\). The set is therefore the epigraph of a convex function and is convex by Definition 1.2.1. The points \((0.25, 4)\) and \((4, 0.25)\) satisfy \(xy = 1 \le 1\), but their midpoint \((2.125, 2.125)\) has \(xy = 4.515625 > 1\), so the set \(\{xy \le 1\}\) fails the chord test. (iii) Both sides are nonnegative, so squaring is allowed, and \(4 + (x - y)^2 \le (x + y)^2\) is \(4 \le 4xy\). The set \(\{(x, y, u) : \lVert (u, x - y) \rVert_2 \le x + y\}\) is a convex cone, the rotated second-order cone, because it is the image of the standard cone \(\{(a, b, c) : \lVert (b, c) \rVert \le a\}\) under an invertible linear map. The slice \(u = 2\) of a convex set is convex. ∎
(One curve, two sides, and what a solver can check) Part (ii) says that one and the same curve, the hyperbola \(xy = 1\), bounds a convex region on one side and a nonconvex region on the other. One side of the curve is convex and the other is not, and the methods for the two sides differ completely. A second-order cone, the set of pairs \((a, b)\) of a scalar and a vector with \(\lVert b \rVert_2 \le a\), is convex, and interior-point methods solve problems over it with the same guarantees as linear programs (Section 4.8). A solver given \(xy \ge 1\) therefore recognizes, by (iii), a conic constraint it can treat exactly, with a certificate. Given \(xy \le 1\) it must relax the product and branch. SCIP 10 lists exactly this recognition as a new feature, with the bullet "extended SOC detection to simple bilinear constraints, e.g., x*y >= 1" in its change log.SCIP, CHANGELOG, release 10.0.0 (24 November 2025), section "Features and Performance Improvements > Nonlinearity", quoted verbatim; github.com/scipopt/scip. The release is described in C. Hojny et al., "The SCIP Optimization Suite 10.0", arXiv 2511.18580 (2025). The rotated second-order cone and what is representable with it are in A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization (SIAM, 2001); Section 4.8 returns to cones. The converse of part (i) is false. A function with convex sublevel sets need not be convex, as \(g(x) = \sqrt{|x|}\) shows, so convexity of the constraint functions is a sufficient condition for convexity of the feasible set and not a characterization. In practice it is the condition that matters, because it is the one a solver can check on the functions it is handed. Even that check has limits. Deciding whether a polynomial of degree four is convex is NP-hard, so solvers accept convexity as declared by the modeller, or detect it by syntactic rules applied to the expression, rather than proving it.A. A. Ahmadi, A. Olshevsky, P. A. Parrilo and J. N. Tsitsiklis, "NP-hardness of deciding convexity of quartic polynomials and related problems", Mathematical Programming 137 (2013). For quadratics the test is positive semidefiniteness of the Hessian and is polynomial. The syntactic rules are the tree walks of R. Fourer, C. Maheshwari, A. Neumaier, D. Orban and H. Schichl, "Convexity and concavity detection in computational graphs: tree walks for convexity assessment", INFORMS Journal on Computing 22 (2010); Section 2.5 describes them.
Proposition 1.2.5: one curve, xy = 1, and its two sides (x, y > 0)
the hyperbola xy = 1
/ \
the side xy >= 1 the side xy <= 1
y >= 1/x, the epigraph of (0.25, 4) and (4, 0.25)
a convex function: convex lie in it; their midpoint
| (2.125, 2.125) has
| (iii): squaring, xy = 4.515625 > 1:
| 4 + (x - y)^2 not convex
| <= (x + y)^2 |
v v
|| (2, x - y) ||_2 <= x + y the solver relaxes the
the slice u = 2 of a rotated product and branches
second-order cone: treated
exactly, with a certificate
(The watershed, and the two kinds of hole) Rockafellar stated the organizing principle in one sentence in 1993: "In fact the great watershed in optimization isn't between linearity and nonlinearity, but convexity and nonconvexity."R. T. Rockafellar, "Lagrange multipliers and optimality", SIAM Review 35 (1993); the sentence is quoted verbatim and without a page number. Mittelmann's recent INFORMS benchmark talks quote it on the slide that introduces the selected benchmarks, which he sorts into a convex group (LP feasibility, MILP, SDP) and a nonconvex group (NLP, MIQCP, MINLP): H. D. Mittelmann, "Latest progress in optimization software", INFORMS Annual Meeting, 28 October 2025, slide 10, plato.asu.edu/talks/informs2025.pdf. Section 5.2 reports those benchmarks. The sentence covers the two kinds of nonconvexity this series is about. Integrality is one kind. The set \(\{0, 1\}\) contains \(0\) and \(1\) but not their midpoint. A feasible set that contains points with \(y = 0\) and points with \(y = 1\) therefore fails the chord test, and so does any feasible set whose points take two or more values of an integer variable. A constraint that bends the wrong way is the other kind, and \(xy \le 1\) is the simplest instance. Section 1.4 shows that the two kinds are one. The condition \(y \in \{0, 1\}\) is the concave quadratic inequality \(y - y^2 \le 0\) on \([0, 1]\), and dropping the integrality is the same operation as replacing that concave term by its chord. A MINLP is a problem with both kinds of hole. Every method of Sections 2 to 4 replaces the nonconvex feasible set by a convex set containing it, solves over the convex set, and then accounts for the difference, by search or by a better formulation.
The figure below draws the two halves of the watershed. On the left is a function of one variable with three valleys, \(f(x) = 0.15x^4 - 1.2x^2 + 0.9\sin 3x + 0.3x\) on \([-3, 3]\), and a gradient descent \(x \leftarrow x - 0.03\, f'(x)\) started from a point the reader can drag. Started at \(x = 0.90\), the descent ends at \(x = 1.68\) with \(f = -2.54\), while the lowest valley sits at \(x = -2.36\) with \(f = -3.38\). The descent cannot tell that it has stopped short. At its stopping point the first-order conditions hold exactly as they do at the global minimizer, and that is the content of the NLP cell of the class table. The basin of a valley is the set of starts whose descent ends there. The control "random starts N" adds \(N\) uniformly drawn starts and colours each by the valley its descent reaches. The figure's statistics compare the share of starts that found the lowest valley with the share of the interval that its basin occupies, which is \(p = 28.8\%\). The other two basins have \(28.5\%\) and \(42.7\%\), and the watersheds between them lie at \(x = -1.27\) and \(x = 0.44\). With \(N = 10\) the figure's seeded starts find the lowest valley \(4\) times, a rate of \(0.40\), and it prints the probability that at least one of \(N\) independent starts does so, \(1 - (1 - p)^N = 0.967\). With \(N = 30\) the count is \(8\) of \(30\) and the probability exceeds \(99.9\%\). With \(N = 100\) it is \(27\) of \(100\), close to \(p\). The guarantee of many starts is probabilistic and never reaches one, which is the point Section 1.3 makes precise about multistart.
On the right is the chord test. One chord runs through the convex disk and one through the set \(xy \le 1\) from \((0.25, 4)\) to \((4, 0.25)\). A slider walks a point along both at once. At \(t = 0.5\) the point on the second chord is \((2.13, 2.13)\) with \(xy = 4.52\), outside the set, while the point on the first chord is inside wherever it is.
The following program recomputes the left panel's numbers with the figure's own rules and checks the chord test of Proposition 1.2.5 on the two sets of the right panel. The rules are a fixed-step descent, the valleys as sign changes of \(f'\) on a grid of 3,000 steps, and the basins by descents from 1,201 equally spaced starts.
# The watershed figure, recomputed.
#
# The left panel: three valleys, a descent that stops in the nearest one,
# and the share of the interval draining into each valley. Then the chord
# test of the right panel.
import numpy as np
f = lambda x: 0.15 * x**4 - 1.2 * x**2 + 0.9 * np.sin(3 * x) + 0.3 * x
df = lambda x: 0.6 * x**3 - 2.4 * x + 2.7 * np.cos(3 * x) + 0.3
L, U = -3.0, 3.0
def descend(x, step=0.03, iters=400):
"""Gradient descent with a fixed step, as the figure runs it."""
for _ in range(iters):
d = df(x)
if abs(d) < 1e-5:
break
x = min(U, max(L, x - step * d))
return x
# valleys: sign changes of f' from - to + on a fine grid;
# watersheds: changes from + to -
g = L + (U - L) * np.arange(3001) / 3000
s = np.sign(df(g))
valleys = [g[i] for i in range(1, 3001) if s[i - 1] < 0 <= s[i]]
sheds = [g[i] for i in range(1, 3001) if s[i - 1] > 0 >= s[i]]
print("valleys (x, f):",
", ".join(f"({v:.3f}, {f(v):.3f})" for v in valleys))
print("watersheds at x =", ", ".join(f"{w:.2f}" for w in sheds))
x_end = descend(0.90)
print(f"descent from x = 0.90 ends at x = {x_end:.2f}, "
f"f = {f(x_end):.2f}; "
f"lowest valley f = {min(f(v) for v in valleys):.2f}")
# basins: a descent from each of 1201 equally spaced starts, counted by
# the valley it reaches
ends = np.array([descend(x)
for x in L + (U - L) * np.arange(1201) / 1200])
share = [np.mean(np.abs(ends - v) < 0.05) for v in valleys]
print("basin shares:", ", ".join(f"{p:.3f}" for p in share))
p = share[int(np.argmin([f(v) for v in valleys]))]
print(f"global basin p = {p:.3f}; 1 - (1 - p)^N for N = 1, 10, 30: "
+ ", ".join(f"{1 - (1 - p)**N:.3f}" for N in (1, 10, 30)))
# the chord test: a point on the chord between two points of x*y <= 1
# leaves the set; a chord of the disk does not
a, b = np.array([0.25, 4.0]), np.array([4.0, 0.25])
m = 0.5 * (a + b)
print(f"chord of x*y <= 1: a = {a}, b = {b}, midpoint {m}, "
f"product {m[0] * m[1]:.6f} > 1: outside")
# two points of the unit disk
c, d = np.array([-0.6, 0.8]), np.array([0.8, 0.6])
t = np.linspace(0, 1, 11)
pts = np.outer(1 - t, c) + np.outer(t, d)
print(f"chord of the disk: max |p| along it "
f"{np.max(np.hypot(pts[:, 0], pts[:, 1])):.3f} <= 1: "
"inside at every t")
It prints the valleys \((-2.356, -3.382)\), \((-0.782, -1.555)\) and \((1.682, -2.540)\), the watersheds \(-1.27\) and \(0.44\), and the descent from \(0.90\) ending at \(1.68\) with \(f = -2.54\). It then prints the basin shares \(0.288\), \(0.285\) and \(0.427\), the probabilities \(0.288\), \(0.967\) and \(1.000\) for \(N = 1, 10, 30\), the midpoint \((2.125, 2.125)\) with product \(4.515625\), and a chord of the disk that stays inside. The cost is one function evaluation per descent step and at most 400 steps per start. The 1,201 descents are independent of one another, so a batch of starts is the simplest parallel pattern in the series: one thread per start, with no communication until the final count. Its weakness is that no number of starts certifies the answer, and that is why the rest of the series is about bounds rather than about starting points.
Where this is used
A solver decides at the outset which side of the watershed each function is on, by the declared type of the problem or by a convexity check on its expression tree. That decision selects the whole algorithm: tangent cuts and local solves with certificates on the convex side, envelopes and spatial branching on the nonconvex side.
What parallelizes
The same decision selects what a GPU is asked to do. In the class table, the LP cell and the convex half of the NLP cell are the ones with certificates, and they are what the bounding step of every method in Sections 6 and 7 batches across many nodes at once. The MILP cell's relaxation is an LP, batched the same way, and its tree is the exponential search of the first difficulty. The nonconvex half of the NLP cell and the MINLP cell are the ones where a parallel machine can only speed up a search that remains, in the worst case, exponential. Everything that is batched across nodes in Sections 6 and 7 is therefore a convex relaxation, because a convex relaxation is the only object whose solution a device can be trusted to return as a bound.
The watershed, decided once per function at the outset
a function of the model
|
the declared type of the problem, or a
convexity check on its expression tree
/ \
convex nonconvex
| |
tangent cuts, local solves envelopes, spatial
with certificates branching
| |
LP, convex half of NLP: nonconvex half of NLP, MINLP:
bounds batched across parallelism only speeds up a
many nodes (Sections 6, 7) worst-case exponential search
What a local method proves
The table of Section 1.1 said that a local method applied to a nonlinear program finds a point satisfying the first-order conditions, and that this is a certificate of local optimality at best. This subsection says what those conditions are and why a minimizer must satisfy them. It then says when they are also sufficient, how the two standard local methods compute a point that satisfies them, and what a run of such a method does and does not prove. It is needed now for two reasons. Every incumbent of a nonlinear problem in this series comes from a local method or is checked by one, so what an incumbent certifies is what a local method certifies. And every bound in this series rests on the Lagrangian defined here, which Section 2.2 turns into the dual function.
Throughout this subsection the problem is continuous,
\[\min_{x \in \mathbb{R}^n}\ f(x) \quad \text{subject to} \quad g_i(x) \le 0, \qquad i = 1, \dots, m, \tag{1.3.1}\]with \(f\) and the \(g_i\) continuously differentiable. Equality constraints \(h_j(x) = 0\) are handled in the same way with multipliers of unrestricted sign, and they are omitted to keep the displays short. Where they appear they are kept as equalities, rather than written as the pair \(h \le 0\), \(-h \le 0\) of Section 1.1, because that pair makes the active gradients dependent and so fails the constraint qualifications LICQ and MFCQ of Definition 1.3.3 below at every feasible point. The qualifications below are stated for the equality form in Nocedal and Wright, Chapter 12, cited at Definition 1.3.3. The feasible set is \(\mathcal F = \{x : g(x) \le 0\}\). Local and global minimizers are as in Definition 1.2.2: a feasible \(\bar x\) is a local minimizer if \(f(\bar x) \le f(x)\) for every feasible \(x\) in some neighbourhood of \(\bar x\), a strict local minimizer if \(f(\bar x) < f(x)\) for every feasible \(x \ne \bar x\) in some such neighbourhood, and a global minimizer if the weak inequality holds for every feasible \(x\).
First-order conditions
Definition 1.3.1 (active set, Lagrangian, KKT point). The active set at a feasible point \(\bar x\) is \(A(\bar x) = \{ i : g_i(\bar x) = 0 \}\), the set of constraints that hold with equality there. The Lagrangian of (1.3.1) is
\[L(x, \lambda) \;=\; f(x) + \sum_{i=1}^m \lambda_i\, g_i(x) \;=\; f(x) + \lambda^\top g(x), \qquad \lambda \in \mathbb{R}^m_{\ge 0},\]and the components of \(\lambda\) are the multipliers. A pair \((\bar x, \bar\lambda)\) is a Karush–Kuhn–Tucker point, or KKT point, of (1.3.1) if
\[\nabla f(\bar x) + \sum_{i=1}^m \bar\lambda_i \nabla g_i(\bar x) = 0, \qquad g(\bar x) \le 0, \qquad \bar\lambda \ge 0, \qquad \bar\lambda_i\, g_i(\bar x) = 0 \quad (i = 1, \dots, m). \tag{1.3.2}\]The four conditions are called stationarity, primal feasibility, dual feasibility and complementary slackness. Complementary slackness says that an inactive constraint carries a zero multiplier, so the stationarity equation involves only the gradients of the active constraints.
Example 1.3.2 (the four conditions on one instance). Consider \(\min (x - 2)^2\) subject to \(x \le 1\), so that \(f(x) = (x - 2)^2\), \(g_1(x) = x - 1\) and \(m = 1\). The Lagrangian is \(L(x, \lambda) = (x - 2)^2 + \lambda (x - 1)\). The four conditions read \(2(x - 2) + \lambda = 0\), \(x \le 1\), \(\lambda \ge 0\) and \(\lambda (x - 1) = 0\). If \(\lambda = 0\), stationarity gives \(x = 2\), which violates \(x \le 1\). So the constraint is active, \(A(\bar x) = \{1\}\), and stationarity at \(x = 1\) gives \(\lambda = 2\). The unique KKT point is \((\bar x, \bar\lambda) = (1, 2)\). Whether it is a minimizer, and what the number \(2\) means, are settled below by Theorem 1.3.4 and Proposition 1.3.5.
Definition 1.3.3 (constraint qualifications). Let \(\bar x\) be feasible. The linear independence constraint qualification (LICQ) holds at \(\bar x\) if the active gradients \(\{\nabla g_i(\bar x) : i \in A(\bar x)\}\) are linearly independent. The Mangasarian–Fromovitz constraint qualification (MFCQ) holds if there is a direction \(d\) with \(\nabla g_i(\bar x)^\top d < 0\) for every \(i \in A(\bar x)\). For a problem whose \(g_i\) are convex, Slater's condition holds if some point \(\tilde x\) has \(g_i(\tilde x) < 0\) for every nonlinear \(g_i\), with \(g_i(\tilde x) \le 0\) allowed for the affine ones. LICQ implies MFCQ, and under MFCQ the multipliers at a local minimizer form a nonempty bounded set.O. L. Mangasarian and S. Fromovitz, "The Fritz John necessary optimality conditions in the presence of equality and inequality constraints", Journal of Mathematical Analysis and Applications 17 (1967); M. Slater, "Lagrange multipliers revisited", Cowles Commission Discussion Paper, Mathematics 403 (1950), reprinted in Traces and Emergence of Nonlinear Programming (Birkhäuser, 2014); J. Gauvin, "A necessary and sufficient regularity condition to have bounded multipliers in nonconvex programming", Mathematical Programming 12 (1977), 136–138, for the boundedness statement; J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer, 2006), Chapter 12, for the relations between the qualifications. A constraint qualification says that the linearized constraints describe the feasible set near \(\bar x\) faithfully. The cusp example after the next theorem shows what happens when they do not.
Theorem 1.3.4 (first-order necessary conditions; Karush 1939; Kuhn and Tucker 1951; John 1948). Let \(f, g_i \in C^1\) and let \(\bar x\) be a local minimizer of (1.3.1). (i) There exist \((\lambda_0, \lambda) \ne 0\) with \(\lambda_0 \ge 0\), \(\lambda \ge 0\), \(\lambda_0 \nabla f(\bar x) + \sum_i \lambda_i \nabla g_i(\bar x) = 0\) and \(\lambda_i g_i(\bar x) = 0\) for all \(i\). (ii) If LICQ or MFCQ holds at \(\bar x\), or if the problem is convex and Slater's condition holds, then \(\lambda_0\) can be taken equal to \(1\): there is a \(\bar\lambda \ge 0\) such that \((\bar x, \bar\lambda)\) is a KKT point.W. Karush, Minima of Functions of Several Variables with Inequalities as Side Conditions, MSc thesis, University of Chicago (1939), reprinted in Traces and Emergence of Nonlinear Programming (Birkhäuser, 2014); H. W. Kuhn and A. W. Tucker, "Nonlinear programming", Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability (University of California Press, 1951), 481–492; F. John, "Extremum problems with inequalities as subsidiary conditions", in Studies and Essays Presented to R. Courant on his 60th Birthday (Interscience, 1948). Full proofs of both parts: Nocedal and Wright (2006), Chapter 12.
Proof sketch of (ii) under LICQ. Let \(T(\bar x)\) be the tangent cone of \(\mathcal F\) at \(\bar x\), the set of limits of directions \((x_k - \bar x)/t_k\) with \(x_k \in \mathcal F\), \(x_k \to \bar x\) and \(t_k \downarrow 0\), and let \(T_{\mathrm{lin}}(\bar x) = \{ d : \nabla g_i(\bar x)^\top d \le 0,\ i \in A(\bar x) \}\) be the linearized cone. Local optimality gives \(\nabla f(\bar x)^\top d \ge 0\) for every \(d \in T(\bar x)\), since a feasible sequence approaching \(\bar x\) along a direction with \(\nabla f(\bar x)^\top d < 0\) would eventually have \(f(x_k) < f(\bar x)\). The inclusion \(T(\bar x) \subseteq T_{\mathrm{lin}}(\bar x)\) always holds. Under LICQ the implicit function theorem gives the reverse inclusion. For every \(d\) in the linearized cone the ray \(\bar x + t d\) can be bent back onto the active constraints by a correction of order \(t^2\), which produces a feasible curve with tangent \(d\). Hence \(\nabla f(\bar x)^\top d \ge 0\) whenever \(\nabla g_i(\bar x)^\top d \le 0\) for all active \(i\). Farkas' lemma says that for a matrix \(G\) and a vector \(c\) exactly one of two systems has a solution: \(G d \le 0\) with \(c^\top d > 0\), or \(G^\top \lambda = c\) with \(\lambda \ge 0\). Applied with the rows \(\nabla g_i(\bar x)^\top\), \(i \in A(\bar x)\), as \(G\) and \(c = -\nabla f(\bar x)\), it turns the implication just proved into the existence of \(\bar\lambda_i \ge 0\), \(i \in A(\bar x)\), with \(-\nabla f(\bar x) = \sum_{i \in A(\bar x)} \bar\lambda_i \nabla g_i(\bar x)\). Setting \(\bar\lambda_i = 0\) off the active set gives (1.3.2). ∎
(What the conditions say geometrically) The picture is this. At a constrained minimizer the direction of steepest descent, \(-\nabla f(\bar x)\), points into the cone spanned by the outward normals of the active constraints, so every direction that decreases \(f\) leaves the feasible set. A constraint qualification guarantees that the gradients of the constraints describe the geometry of \(\mathcal F\) near \(\bar x\) faithfully. Without one the conditions can degenerate to \(\lambda_0 = 0\), and then they say nothing about \(f\). The standard example is \(\min x_1\) subject to \(x_2 \le x_1^3\) and \(x_2 \ge 0\). The origin is the minimizer. The two active gradients there are \((0, 1)\) and \((0, -1)\), which are dependent, and no \(\lambda \ge 0\) makes \((1, 0) + \lambda_1 (0, 1) + \lambda_2 (0, -1)\) vanish. Only the Fritz John form (i) holds, with \(\lambda_0 = 0\) and \(\lambda_1 = \lambda_2 = 1\). A solver that assumes a constraint qualification stalls at such a point.
(Multipliers as prices) The multipliers have a second meaning, which Rockafellar's paper on Lagrange multipliers places at the centre of the subject: they are the prices of the constraints.R. T. Rockafellar, SIAM Review 35 (1993), cited in Section 1.2. Write \(v(u) = \inf\{ f(x) : g(x) \le u \}\) for the optimal value as a function of the right-hand sides, so that \(v(0) = z^\star\). Loosening constraint \(i\) by \(u_i\) is worth about \(\bar\lambda_i u_i\). For a convex problem this is a theorem that comes with the global optimality of KKT points.
Proposition 1.3.5 (KKT points of convex problems; multipliers as prices). Let \(f\) and every \(g_i\) be convex and differentiable, and let \((\bar x, \bar\lambda)\) be a KKT point of (1.3.1). Then (i) \(\bar x\) is a global minimizer, and (ii) for every \(u \in \mathbb{R}^m\),
\[v(u) \;\ge\; v(0) - \bar\lambda^\top u ,\]so that \(-\bar\lambda\) is a subgradient of \(v\) at \(0\), and \(\partial v / \partial u_i\,(0) = -\bar\lambda_i\) wherever \(v\) is differentiable.
Proof. Let \(x\) satisfy \(g(x) \le u\). Convexity of \(f\), then stationarity, then convexity of each \(g_i\) multiplied by \(\bar\lambda_i \ge 0\), then complementary slackness, then \(g(x) \le u\) give
\[f(x) \;\ge\; f(\bar x) + \nabla f(\bar x)^\top (x - \bar x) \;=\; f(\bar x) - \sum_i \bar\lambda_i \nabla g_i(\bar x)^\top (x - \bar x) \;\ge\; f(\bar x) - \sum_i \bar\lambda_i \big( g_i(x) - g_i(\bar x) \big) \;=\; f(\bar x) - \bar\lambda^\top g(x) \;\ge\; f(\bar x) - \bar\lambda^\top u .\]Taking the infimum over such \(x\) gives \(v(u) \ge f(\bar x) - \bar\lambda^\top u\). With \(u = 0\) this says \(f(x) \ge f(\bar x)\) for every feasible \(x\), which is (i) and also \(v(0) = f(\bar x)\). Substituting back gives (ii). ∎
On Example 1.3.2 the proposition says two things. The problem is convex, so the KKT point \(x = 1\) is the global minimizer, with value \(1\). And the multiplier is the price of the constraint. Loosening it to \(x \le 1 + \delta\) gives the value \(v(\delta) = (1 + \delta - 2)^2 = 1 - 2\delta + \delta^2\), a decrease of \(2\delta - \delta^2 \approx \bar\lambda \delta\), which is \(0.0199\) at \(\delta = 0.01\).
For a nonconvex problem the price interpretation survives locally, as the sensitivity theorem of Fiacco and McCormick. It needs LICQ, strict complementarity, which means that every active constraint carries a strictly positive multiplier, and the second-order sufficient condition of Theorem 1.3.8 below.A. V. Fiacco and G. P. McCormick, Nonlinear Programming: Sequential Unconstrained Minimization Techniques (Wiley, 1968; SIAM Classics reprint, 1990); S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004), Section 5.6, for the convex case stated above. Global optimality does not survive, and the next example shows what is lost.
Example 1.3.6 (a KKT point is a candidate, not a proof). Change one sign in Example 1.3.2 and widen the interval: \(\min -(x - 2)^2\) on \([0, 3]\), with multipliers \(\mu_0\) for \(-x \le 0\) and \(\mu_1\) for \(x - 3 \le 0\). Stationarity reads \(-2(x - 2) - \mu_0 + \mu_1 = 0\), and there are three KKT points. At \(x = 2\) both multipliers vanish and the value is \(0\), which is the maximum of the objective on the interval. At \(x = 0\) the first constraint is active with \(\mu_0 = 4\) and the value is \(-4\). At \(x = 3\) the second constraint is active with \(\mu_1 = 2\) and the value is \(-1\). A descent method started at \(x = 2.5\) moves to the right, because \(-f'(x) = 2(x - 2) > 0\) there, stops at \(x = 3\), verifies the KKT conditions, and reports \(-1\). The report is correct: \(x = 3\) is a strict local minimizer. The answer is wrong by a factor of four, because the global minimizer is \(x = 0\). Nothing computed at \(x = 3\) reveals this. The two problems are quadratics over intervals, one with a single KKT point and one with three, and the same local method handles both. The difference is the sign of one coefficient, which is what Section 1.2 called the watershed.
Second-order conditions, and what can be certified
The KKT conditions are first-order statements. The second-order information a solver has at termination is the curvature of the Lagrangian along the directions that keep the active constraints active.
Definition 1.3.7 (critical cone). At a KKT point \((\bar x, \bar\lambda)\) the critical cone is
\[C(\bar x, \bar\lambda) \;=\; \{ d : \nabla g_i(\bar x)^\top d = 0 \ \text{for } i \in A(\bar x) \text{ with } \bar\lambda_i > 0, \quad \nabla g_i(\bar x)^\top d \le 0 \ \text{for } i \in A(\bar x) \text{ with } \bar\lambda_i = 0 \}.\]Theorem 1.3.8 (second-order conditions; McCormick 1967). Let \(f, g_i \in C^2\), let \((\bar x, \bar\lambda)\) be a KKT point, and write \(\nabla^2_{xx} L = \nabla^2 f(\bar x) + \sum_i \bar\lambda_i \nabla^2 g_i(\bar x)\). (i) Necessary (the second-order necessary condition, SONC): if \(\bar x\) is a local minimizer and LICQ holds at \(\bar x\), then \(d^\top \nabla^2_{xx} L\, d \ge 0\) for every \(d \in C(\bar x, \bar\lambda)\). (ii) Sufficient (the second-order sufficient condition, SOSC): if \(d^\top \nabla^2_{xx} L\, d > 0\) for every nonzero \(d \in C(\bar x, \bar\lambda)\), then \(\bar x\) is a strict local minimizer, and no constraint qualification is needed.G. P. McCormick, "Second order conditions for constrained minima", SIAM Journal on Applied Mathematics 15 (1967); proofs in Nocedal and Wright (2006), Section 12.5.
(Where the curvature lives) It is the Lagrangian and not the objective that carries the curvature, because the constraints bend, and the test is needed only along directions that keep the strongly active constraints active. When strict complementarity holds, so that \(\bar\lambda_i > 0\) for every active \(i\), the critical cone is the null space of the active gradients. The test is then the definiteness of a projected Hessian, which is a matrix computation. When some active multiplier vanishes the cone has flat sides, the test becomes the nonnegativity of a quadratic form on a cone, and that question is hard in its own right.
Theorem 1.3.9 (checking local optimality is hard; Murty and Kabadi 1987; Pardalos and Schnitger 1988; Ahmadi and Zhang 2022). (i) Deciding whether a symmetric integer matrix \(Q\) is copositive, that is, whether \(x^\top Q x \ge 0\) for all \(x \ge 0\), is co-NP-complete. (ii) Deciding whether a given feasible point is a local minimizer of a quadratic program \(\min\{ x^\top H x + c^\top x : Ax \le b \}\) is NP-hard, and this remains true for constrained problems under the restrictions studied by Pardalos and Schnitger. (iii) Unless P = NP, no polynomial-time algorithm finds a point within Euclidean distance \(c^n\) of a local minimizer of an \(n\)-variable quadratic function over a polytope, for any constant \(c \ge 0\).K. G. Murty and S. N. Kabadi, "Some NP-complete problems in quadratic and nonlinear programming", Mathematical Programming 39 (1987); P. M. Pardalos and G. Schnitger, "Checking local optimality in constrained quadratic programming is NP-hard", Operations Research Letters 7 (1988); A. A. Ahmadi and J. Zhang, "On the complexity of finding a local minimizer of a quadratic function over a polytope", Mathematical Programming 195 (2022).
Proof sketch of (ii) from (i). Take \(c = 0\), the constraints \(-x \le 0\) and the point \(\bar x = 0\). The feasible set is the cone \(\mathbb{R}^n_{\ge 0}\) and the objective \(x^\top H x\) is homogeneous of degree two, so \(0\) is a local minimizer if and only if \(x^\top H x \ge 0\) for all \(x \ge 0\), which is the copositivity of \(H\). ∎
(What a local method can certify) The consequence for every solver in this series is that a local method certifies second-order sufficiency only in the nondegenerate case. In the degenerate case it certifies a KKT point and stops. The sentence "local optimization is easy and global optimization is hard" is therefore a statement about typical behaviour, not a theorem. What can be checked independently of the solver, for any candidate point it returns, is the following.
Algorithm 1.3.10 (verifying a candidate point).
Input a point x_bar; f, g_1..g_m with first and second derivatives;
tolerances eps_feas, eps_stat.
Output one of {not feasible; feasible but not KKT; KKT, degenerate
(undecided); KKT, not a local minimizer; undecided (zero
curvature); KKT + SOSC: strict local minimizer}.
1. if max_i g_i(x_bar) > eps_feas:
return "not feasible", with the violation.
2. A := { i : |g_i(x_bar)| <= eps_feas };
G := the matrix whose rows are grad g_i(x_bar), i in A.
3. lambda_A := argmin_{lambda >= 0} || grad f(x_bar) + G^T lambda ||_2
(a nonnegative least-squares problem of size |A|);
lambda_i := 0 for i not in A.
4. if || grad f(x_bar) + G^T lambda_A || > eps_stat:
return "feasible, not KKT", with the stationarity residual.
5. if some lambda_i = 0 for i in A, or rank G < |A|:
return "KKT, degenerate: the second-order test is a
copositivity question (Theorem 1.3.9); report the projected
Hessian spectrum as information only".
6. Z := an orthonormal basis of the null space of G;
H := Z^T [grad^2 f(x_bar) + sum_i lambda_i grad^2 g_i(x_bar)] Z.
7. if lambda_min(H) > 0: return "KKT + SOSC: strict local minimizer";
if lambda_min(H) < 0:
return "KKT, not a local minimizer", with the descent
direction Z v for the eigenvector v;
else return "undecided".
Invariant
every verdict is implied by Theorems 1.3.4 and 1.3.8, so a verdict
never depends on how x_bar was produced.
The cost of one verdict is one evaluation of the first and second derivatives, one nonnegative least-squares problem of size \(|A|\) and one symmetric eigendecomposition of order \(n - |A|\). The verdicts for different candidates are independent and share the same control flow, so a pool of incumbents or multistart results is checked in one batch, and that batch is the natural verification kernel of a GPU solver.
Example 1.3.11 (two KKT points, one of them a saddle). Let \(f(x) = (x_1 - \tfrac12)^2 + x_2^2\) and \(g(x) = 1 - x_1^2 - x_2^2 \le 0\): the point of the closed exterior of the unit disc nearest to \((\tfrac12, 0)\). The feasible set is nonconvex, since it is the complement of an open disc. The unconstrained minimizer \((\tfrac12, 0)\) is infeasible, so the constraint is active at every KKT point, and on the circle \(x = (\cos t, \sin t)\) stationarity \(\nabla f \parallel \nabla g\) reduces to \(\sin t = 0\). There are exactly two KKT points: \((1, 0)\) with \(\lambda = \tfrac12\) and \((-1, 0)\) with \(\lambda = \tfrac32\). The Hessian of the Lagrangian is \((2 - 2\lambda) I\). At \((1, 0)\) it is \(+I\), the second-order sufficient condition holds, and the point is the global minimizer with \(f = \tfrac14\). At \((-1, 0)\) it is \(-I\) and the curvature along the tangent direction \((0, 1)\) is \(-1\), so the necessary condition fails. The point is a saddle of the constrained problem with \(f = \tfrac94\): moving along the circle decreases \(f\) and moving outward increases it. The program below runs Newton's method on the KKT system (1.3.2), with the constraint held active, from three starting points, and applies the second-order test of Algorithm 1.3.10 to whatever it converges to. It then prints the arithmetic of Examples 1.3.2 and 1.3.6.
# What a local method proves.
#
# The KKT points of Example 1.3.11, the nearest point to (1/2, 0) outside
# the unit disc:
# min (x1 - 1/2)^2 + x2^2 s.t. g(x) = 1 - x1^2 - x2^2 <= 0.
# There are two KKT points; Newton on the KKT system reaches whichever is
# nearer, and the second-order test tells them apart. Then the
# one-variable pair of Examples 1.3.2 and 1.3.6.
import numpy as np
np.set_printoptions(precision=4, suppress=True)
gf = lambda x: np.array([2 * (x[0] - 0.5), 2 * x[1]]) # grad f
gg = lambda x: np.array([-2 * x[0], -2 * x[1]]) # grad g
f = lambda x: (x[0] - 0.5)**2 + x[1]**2
g = lambda x: 1 - x[0]**2 - x[1]**2
def newton_kkt(x, lam, iters=50):
"""Newton on F(x, lam) = [grad f + lam grad g ; g] = 0.
The constraint is held active.
"""
for k in range(iters):
F = np.concatenate([gf(x) + lam * gg(x), [g(x)]])
if np.linalg.norm(F) < 1e-12:
return x, lam, k
J = np.block([[(2 - 2 * lam) * np.eye(2), gg(x)[:, None]],
[gg(x)[None, :], np.zeros((1, 1))]])
# least squares guards against a singular Jacobian
step = np.linalg.lstsq(J, -F, rcond=None)[0]
x, lam = x + step[:2], lam + step[2]
return x, lam, iters
for x0 in [(0.9, 0.1), (-0.9, 0.1), (-0.6, 0.8)]:
x, lam, k = newton_kkt(np.array(x0, float), 1.0)
# print -0.0 as 0
x = np.where(np.abs(x) < 1e-9, 0.0, x)
# the tangent to the circle: the critical cone
d = np.array([-x[1], x[0]])
# d' [Hess f + lam Hess g] d
curv = d @ ((2 - 2 * lam) * np.eye(2)) @ d
if curv > 0:
verdict = "strict local minimizer (SOSC holds)"
else:
verdict = "NOT a local minimizer (SONC fails)"
print(f"start {x0}: {k:2d} Newton steps -> x = {x}, "
f"lambda = {lam:.4f},")
print(f" f = {f(x):.4f}, "
f"stationarity {np.linalg.norm(gf(x) + lam * gg(x)):.1e},")
print(f" tangent curvature {curv:+.2f}: {verdict}")
# Examples 1.3.2 and 1.3.6: min (x-2)^2 s.t. x <= 1, and its nonconvex
# twin min -(x-2)^2 on [0, 3]
delta = 0.01
print()
print("convex: optimum x = 1, lambda = 2;")
print(f" value at x <= 1 + delta: {(1 + delta - 2)**2:.4f} "
f"vs {1.0:.4f},")
print(f" drop {1 - (1 + delta - 2)**2:.4f} = 2 delta - delta^2")
h = lambda x: -(x - 2)**2
kkt = {0.0: ("mu_0 = 4", h(0.0)),
2.0: ("mu = 0", h(2.0)),
3.0: ("mu_1 = 2", h(3.0))}
for x, (m, v) in kkt.items():
print(f"nonconvex twin: KKT point x = {x:.0f} ({m}), "
f"value {v + 0.0:+.0f}")
# projected gradient descent from 2.5
x = 2.5
for _ in range(200):
x = min(3.0, max(0.0, x - 0.1 * (-2 * (x - 2))))
print(f"descent from 2.5 stops at x = {x:.0f} with value {h(x):+.0f};")
print(f" the global minimum is "
f"{min(v for _, v in kkt.values()):+.0f} at x = 0")
Output:
start (0.9, 0.1): 5 Newton steps -> x = [1. 0.], lambda = 0.5000,
f = 0.2500, stationarity 2.2e-16,
tangent curvature +1.00: strict local minimizer (SOSC holds)
start (-0.9, 0.1): 5 Newton steps -> x = [-1. 0.], lambda = 1.5000,
f = 2.2500, stationarity 0.0e+00,
tangent curvature -1.00: NOT a local minimizer (SONC fails)
start (-0.6, 0.8): 7 Newton steps -> x = [-1. 0.], lambda = 1.5000,
f = 2.2500, stationarity 0.0e+00,
tangent curvature -1.00: NOT a local minimizer (SONC fails)
convex: optimum x = 1, lambda = 2;
value at x <= 1 + delta: 0.9801 vs 1.0000,
drop 0.0199 = 2 delta - delta^2
nonconvex twin: KKT point x = 0 (mu_0 = 4), value -4
nonconvex twin: KKT point x = 2 (mu = 0), value +0
nonconvex twin: KKT point x = 3 (mu_1 = 2), value -1
descent from 2.5 stops at x = 3 with value -1;
the global minimum is -4 at x = 0
Newton's method on the KKT system converges to whichever KKT point is nearer. From \((-0.9, 0.1)\) and from \((-0.6, 0.8)\) it converges, with stationarity residual \(0\) to machine precision, to the saddle \((-1, 0)\), and a solver that stopped on the first-order residual alone would report \(2.25\) as its answer. The second-order test, one \(2 \times 2\) curvature, is what separates the two verdicts. The cost here is a \(3 \times 3\) linear solve per iteration. In general it is the factorization discussed next. The three runs are independent, and on a GPU they would be one batch of three.
How a local method computes a KKT point
Three families of local method compute KKT points, and a solver's choice among them is a choice of linear algebra. What the series needs from each is its cost per step, a factorization or a few products; how it behaves when restarted on a nearby problem, which decides its usefulness at the nodes of a tree; and which of its kernels a parallel machine can take over. The first two families are described in some detail because they are what every solver of Section 5 runs, and the third is named because it is the one that fits a GPU.
Sequential quadratic programming applies Newton's method to the KKT system (1.3.2) while keeping track of which constraints are active. At the iterate \((x_k, \lambda_k)\) it solves the quadratic program
\[\min_d\ \nabla f(x_k)^\top d + \tfrac12\, d^\top \nabla^2_{xx} L(x_k, \lambda_k)\, d \quad \text{subject to} \quad g_i(x_k) + \nabla g_i(x_k)^\top d \le 0, \qquad i = 1, \dots, m, \tag{1.3.3}\]whose constraints are the linearized constraints and whose quadratic term is the curvature of the Lagrangian, and sets \(x_{k+1} = x_k + d\) with \(\lambda_{k+1}\) the multipliers of the quadratic program. Near a KKT point at which LICQ, strict complementarity and the second-order sufficient condition hold, the quadratic program identifies the correct active set in finitely many steps. From then on the iteration is Newton's method on the equations of the active set, with quadratic convergence (Robinson 1974). With a quasi-Newton approximation of \(\nabla^2_{xx} L\) the rate is superlinear (Powell 1978).S. M. Robinson, "Perturbed Kuhn–Tucker points and rates of convergence for a class of nonlinear-programming algorithms", Mathematical Programming 7 (1974); M. J. D. Powell, "A fast algorithm for nonlinearly constrained optimization calculations", in Numerical Analysis (Dundee 1977), Lecture Notes in Mathematics 630 (Springer, 1978); S. P. Han, "A globally convergent method for nonlinear programming", Journal of Optimization Theory and Applications 22 (1977), for the line search on the \(\ell_1\) penalty function; R. Fletcher and S. Leyffer, "Nonlinear programming without a penalty function", Mathematical Programming 91 (2002), for the filter; P. T. Boggs and J. W. Tolle, "Sequential quadratic programming", Acta Numerica 4 (1995), for a survey. filterSQP and SNOPT (P. E. Gill, W. Murray and M. A. Saunders, "SNOPT: an SQP algorithm for large-scale constrained optimization", SIAM Review 47 (2005)) are the reference implementations. Far from a solution the step is accepted only if it improves a merit function, such as \(f(x) + \sigma \sum_i \max(0, g_i(x))\) (Han 1977), or if it improves either the objective or the constraint violation, which is Fletcher and Leyffer's filter. Under standard assumptions the iterates converge to KKT points, and the limit is a KKT point and nothing more. The cost of a step is the solution of the quadratic program (1.3.3), itself a sequence of factorizations of the KKT matrix of the current active set. SQP warm-starts well: the next solve can begin from the previous active set and the factors of its KKT matrix, both of which carry over to a nearby problem. That is why it is the natural local solver at the nodes of a tree.
One SQP run: the loop, its two phases, the warm start at a node
+--> (x_k, lambda_k)
| |
| | solve the QP (1.3.3): the linearized constraints,
| | the curvature of the Lagrangian as quadratic term
| v
| a step d and the QP's multipliers
| |
| | accept d
| v
+--- x_{k+1} = x_k + d, lambda_{k+1} = the QP's multipliers
far from a solution near a KKT point with LICQ,
strict complementarity and SOSC
o-------o-------o-------o--------o-----o---o--o-o-> a KKT point,
nothing more
d accepted only if it the QP identifies the correct
improves a merit function active set in finitely many
(Han), or the objective or steps; then Newton on its
the constraint violation equations: quadratic, or
(Fletcher and Leyffer's superlinear with quasi-Newton
filter)
o an iterate x_k; at a nearby problem, such as the next node of
a tree, the solve warm-starts from the previous active set
and the factors of its KKT matrix
Interior-point methods replace the combinatorial question of which constraints are active by a continuous homotopy. Write the problem with slacks as \(\min f(x)\) subject to \(g(x) + s = 0\), \(s \ge 0\), and consider the barrier subproblem \(\min f(x) - \mu \sum_i \ln s_i\) subject to \(g(x) + s = 0\) for a parameter \(\mu > 0\). Its first-order conditions are
\[\nabla f(x) + \nabla g(x)^\top \lambda = 0, \qquad g(x) + s = 0, \qquad S \lambda = \mu e, \qquad s, \lambda > 0, \tag{1.3.4}\]with \(S = \operatorname{diag}(s)\) and \(e\) the vector of ones: the KKT conditions with complementary slackness relaxed from \(s_i \lambda_i = 0\) to \(s_i \lambda_i = \mu\). As \(\mu \downarrow 0\) the solutions \((x(\mu), s(\mu), \lambda(\mu))\) trace the central path to a KKT point (Fiacco and McCormick 1968).
On Example 1.3.2 every quantity of the central path is explicit. The barrier subproblem is \(\min (x - 2)^2 - \mu \ln(1 - x)\), its minimizer is \(x(\mu) = (3 - \sqrt{1 + 2\mu})/2\), the slack is \(s = 1 - x(\mu)\) and the multiplier is \(\lambda = \mu / s = 1 + \sqrt{1 + 2\mu}\). The pair \((x(\mu), \lambda(\mu))\) satisfies the stationarity equation \(2(x - 2) + \lambda = 0\) exactly and complementary slackness only to \(s\lambda = \mu\), and as \(\mu\) falls it slides along the stationarity line to the KKT point \((1, 2)\).
(The Newton system of an interior-point step) Each step is a Newton step on (1.3.4) in the unknowns \((\Delta x, \Delta s, \Delta \lambda)\). After eliminating \(\Delta s\) the Newton system is the symmetric indefinite system
\[\begin{pmatrix} W + \delta_w I & \nabla g(x)^\top \\ \nabla g(x) & -\Sigma^{-1} \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta \lambda \end{pmatrix} \;=\; - \begin{pmatrix} \nabla f(x) + \nabla g(x)^\top \lambda \\ g(x) + s - \Lambda^{-1} (S\lambda - \mu e) \end{pmatrix}, \qquad W = \nabla^2 f(x) + \sum_i \lambda_i \nabla^2 g_i(x), \quad \Lambda = \operatorname{diag}(\lambda), \quad \Sigma = S^{-1} \Lambda ,\]of order \(n + m\), followed by \(\Delta s = -(g(x) + s) - \nabla g(x)\, \Delta x\). The right-hand side is built from the residuals of (1.3.4), and \(\delta_w \ge 0\) is a regularization parameter, zero unless the test that follows fails. The step is a descent direction for the barrier problem only if the matrix has exactly \(n\) positive and \(m\) negative eigenvalues. The inertia of a symmetric matrix is its triple of counts of positive, negative and zero eigenvalues, so the requirement is the inertia \((n, m, 0)\). That inertia is the one a saddle point of the Lagrangian has, a minimum in \(x\) and a maximum in \(\lambda\). The reason is that the lower right block is negative definite, so the matrix has inertia \((n, m, 0)\) exactly when its Schur complement, the matrix \(W + \delta_w I + \nabla g(x)^\top \Sigma\, \nabla g(x)\) that remains after \(\Delta \lambda\) is eliminated, is positive definite. On the central path, where \(\lambda_i s_i = \mu\), that Schur complement with \(\delta_w = 0\) is the Hessian of the barrier function \(f(x) - \mu \sum_i \ln(-g_i(x))\), and a Newton step taken with a positive definite Hessian is a descent direction for it.
(The factorization is the sequential kernel) The solver reads the inertia off the factorization and increases \(\delta_w\) when it is wrong. This is why the inner linear algebra of IPOPT, KNITRO and their relatives is a sparse symmetric indefinite \(LDL^\top\) factorization with dynamic pivoting (MA27, MA57, MUMPS, Pardiso). It is also why that factorization is the sequential kernel of nonlinear optimization. The pivot order is decided during the factorization for stability. The dense fronts, which are the blocks of dense arithmetic inside the factorization, are small, and the arithmetic intensity is low (both terms are defined in Section 7.1).A. Wächter and L. T. Biegler, "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming", Mathematical Programming 106 (2006), which proves convergence to a KKT point of the barrier problem, or to a stationary point of the constraint violation when the problem is locally infeasible; R. H. Byrd, J. Nocedal and R. A. Waltz, "KNITRO: an integrated package for nonlinear optimization", in Large-Scale Nonlinear Optimization (Springer, 2006); S. Mehrotra, "On the implementation of a primal-dual interior point method", SIAM Journal on Optimization 2 (1992), for the predictor–corrector variant that linear programming codes use; I. S. Duff, "MA57: a code for the solution of sparse symmetric definite and indefinite systems", ACM Transactions on Mathematical Software 30 (2004).
(The linear case, and the augmented Lagrangian family) In the linear case \(W = 0\), eliminating one more block gives the normal equations, a system of the form \(A D A^\top \Delta y = r\). Here \(A\) is the constraint matrix, \(D\) is a positive diagonal matrix built from \(s\) and \(\lambda\), and \(r\) is the corresponding residual. That system is positive definite, its sparsity pattern is fixed for the whole run, and a Cholesky factorization without pivoting handles it. This is why barrier methods for linear programming port to GPUs and general nonlinear barrier methods do not without a change of formulation. Section 7.5 returns to that change, the condensed KKT systems of Pacaud, Shin and their co-authors, and to what it costs in conditioning and accuracy.F. Pacaud, S. Shin, A. Montoison, M. Schanen and M. Anitescu, "Condensed interior-point methods for scalable nonlinear programming on GPUs", Mathematical Programming Computation (2026); arXiv 2405.14236. S. Shin, M. Anitescu and F. Pacaud, "Accelerating optimal power flow with GPUs: SIMD abstraction of nonlinear programs and condensed-space interior-point methods", Electric Power Systems Research 236 (2024), 110651. A third family, the augmented Lagrangian methods of LANCELOT and ALGENCAN, minimizes \(f(x) + \tfrac{\rho}{2} \sum_i \max(0, g_i(x) + \lambda_i/\rho)^2\) over the box by a gradient method and updates \(\lambda\) by a first-order step. It needs no factorization at all, only products and projections, at the price of a linear rather than quadratic rate.A. R. Conn, N. I. M. Gould and Ph. L. Toint, "A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds", SIAM Journal on Numerical Analysis 28 (1991); R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, "On augmented Lagrangian methods with general lower-level constraints", SIAM Journal on Optimization 18 (2008).
An interior-point step, from (1.3.4) down to the factorization
a Newton step on (1.3.4) in (dx, ds, dlambda)
|
| eliminate ds; afterwards ds = -(g(x) + s) - grad g(x) dx
v
[ W + delta_w I grad g(x)^T ] [ dx ] a right-hand side
[ grad g(x) -Sigma^-1 ] [ dlambda ] = built from the
residuals
symmetric indefinite, of order n + m
| |
| general W | the linear case, W = 0:
| | eliminate one more block
v v
LDL^T with dynamic pivoting A D A^T dy = r: positive
(MA27, MA57, MUMPS, Pardiso); definite, sparsity fixed for
the step descends only if the the whole run; Cholesky without
inertia read off it is pivoting: why barrier methods
(n, m, 0), else delta_w grows: for LP port to GPUs
the sequential kernel; no GPU
port without a change of
formulation
(What the termination test proves) Whichever of the three a solver runs, its termination test is the same: the three residuals of (1.3.2), stationarity, feasibility and complementarity, below a tolerance such as \(10^{-8}\). The output is an \(\varepsilon\)-KKT point. Its second-order status is available only implicitly, through the inertia of the last factorization in an interior-point method or the active-set Hessian of the last quadratic program in SQP, and its global status is not available at all.
Multistart
The one thing a local method can do about global optimality is to run many times. Let the local method be deterministic, so that every start \(x_0\) in a compact set \(S\) ends at one of finitely many local minimizers, and let the basin of a minimizer be as in Section 1.2. Write \(p\) for the fraction of the volume of \(S\) occupied by the basin of the global minimizer, which is the figure's \(p\).
Proposition 1.3.12 (pure multistart and the Bayesian estimate). (i) With \(N\) independent uniform starts, the probability that none of them lands in the basin of the global minimizer is \((1 - p)^N\). It tends to zero as \(N\) grows and is never zero: multistart has no finite certificate. (ii) Suppose a uniform prior on the number \(W\) of local minimizers and a uniform (Dirichlet) prior on the basin volumes. After \(N\) local searches have found \(w\) distinct minimizers, the posterior expectation of \(W\) is
\[\mathbb{E}[W \mid w, N] \;=\; \frac{w\,(N - 1)}{N - w - 2} \qquad (N > w + 2),\]and the posterior expectation of the total volume of the basins already found is \(\dfrac{(N - w - 1)(N + w)}{N (N - 1)}\). The usual stopping rule is to stop when the first estimate is below \(w + \tfrac12\).Part (ii) is due to C. G. E. Boender and A. H. G. Rinnooy Kan, "Bayesian stopping rules for multistart global optimization methods", Mathematical Programming 37 (1987), and C. G. E. Boender, A. H. G. Rinnooy Kan, G. T. Timmer and L. Stougie, "A stochastic method for global optimization", Mathematical Programming 22 (1982); the formulas and the thresholds \(0.5\) and \(0.995\) are printed in R. H. Byrd, C. L. Dert, A. H. G. Rinnooy Kan and R. B. Schnabel, "Concurrent stochastic methods for global optimization", Mathematical Programming 46 (1990), pp. 6–7. The clustering methods that start one local search per basin are A. H. G. Rinnooy Kan and G. T. Timmer, "Stochastic global optimization methods. Part I: Clustering methods" and "Part II: Multi level methods", Mathematical Programming 39 (1987).
Part (i) is the multiplication rule for independent events. The point of part (ii) is that "how many minimizers have I seen" can be converted into "how many are there", with a prior. Take a function with three basins occupying \(20\%\), \(50\%\) and \(30\%\) of the interval, the deepest valley first. One start finds it with probability \(0.2\), ten starts miss it with probability \(0.8^{10} = 0.107\), and thirty starts miss it with probability \(0.8^{30} = 0.0012\). The watershed figure of Section 1.2 runs this experiment on its three-valley function with \(p = 0.288\) and prints \(1 - (1 - p)^N\) for the \(N\) it is given. Multistart is Monte Carlo integration of the indicator of the global basin. Its guarantee is probabilistic, its cost is one local solve per start, and the solves are independent, which is the first embarrassingly parallel computation in this series. The clustering refinements of Rinnooy Kan and Timmer start a local search from a sample point only if no better sample lies within a critical distance of it, so that in the limit one local search is spent per basin. The sampling and the local solves remain independent, and only the global count \(w\) and the nearest-neighbour pass are shared.
The verification of Algorithm 1.3.10 is where a batch of candidates, from a multistart or from a pool of incumbents, is turned into verdicts. The following C++23 program does it for Example 1.3.11 over a pool of 2,005 candidate points, five chosen by hand and two thousand spread around the circle, with \(T\) threads, each taking every \(T\)-th candidate of the pool. A GPU version has one thread per candidate and the same body.
// Verifying what a local solver returned, for a batch of candidates.
//
// The GPU-shaped kernel of Section 1.3. The problem is
// min f(x) = (x1 - 1/2)^2 + x2^2 s.t. g(x) = 1 - x1^2 - x2^2 <= 0.
// For each candidate point: feasibility, the least-squares multiplier,
// the stationarity residual, and the curvature of the Lagrangian along
// the tangent of the active constraint (the second-order test).
#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdio>
#include <thread>
#include <vector>
// code: 0 infeasible, 1 not KKT, 2 KKT saddle, 3 strict local min
struct Verdict {
int code;
double lambda, resid, curv;
};
static Verdict classify(double x1, double x2, double eps_feas = 1e-8,
double eps_stat = 1e-6) {
const double g = 1.0 - x1 * x1 - x2 * x2;
if (g > eps_feas)
return {0, 0.0, 0.0, 0.0};
const std::array<double, 2> gf{2.0 * (x1 - 0.5), 2.0 * x2};
const std::array<double, 2> gg{-2.0 * x1, -2.0 * x2};
double lambda = 0.0;
if (std::fabs(g) <= eps_feas) {
// active: lambda = argmin_{l >= 0} |gf + l gg|
const double num = -(gf[0] * gg[0] + gf[1] * gg[1]);
const double den = gg[0] * gg[0] + gg[1] * gg[1];
lambda = std::fmax(0.0, num / den);
}
const double r0 = gf[0] + lambda * gg[0];
const double r1 = gf[1] + lambda * gg[1];
const double resid = std::hypot(r0, r1);
if (resid > eps_stat)
return {1, lambda, resid, 0.0};
// interior stationary point: Hess f = 2I > 0
if (std::fabs(g) > eps_feas)
return {3, 0.0, resid, 2.0};
// d^T (Hess f + lambda Hess g) d, |d| = 1
const double curv = 2.0 - 2.0 * lambda;
return {curv > 1e-9 ? 3 : 2, lambda, resid, curv};
}
int main() {
std::vector<std::array<double, 2>> cand{
{1.0, 0.0}, {-1.0, 0.0}, {0.6, 0.8}, {0.5, 0.0}, {2.0, 0.0}};
// a pool of points on the circle
for (int k = 0; k < 2000; ++k) {
const double t = 6.283185307179586 * (k + 0.5) / 2000.0;
cand.push_back({std::cos(t), std::sin(t)});
}
std::vector<Verdict> out(cand.size());
std::array<std::atomic<int>, 4> counts{};
{
const unsigned T =
std::max(1u, std::thread::hardware_concurrency());
// every T-th candidate per thread; on a GPU, one thread per candidate
std::vector<std::jthread> pool;
for (unsigned t = 0; t < T; ++t)
pool.emplace_back([&, t] {
for (std::size_t i = t; i < cand.size(); i += T) {
out[i] = classify(cand[i][0], cand[i][1]);
counts[out[i].code].fetch_add(
1, std::memory_order_relaxed);
}
});
} // jthreads join here
const char* name[] = {"infeasible", "feasible, not KKT",
"KKT, not a local minimizer",
"KKT + SOSC: strict local minimizer"};
for (std::size_t i = 0; i < 5; ++i) {
std::printf("(%+.2f, %+.2f): %s\n",
cand[i][0], cand[i][1], name[out[i].code]);
std::printf(" lambda %.3f residual %.1e"
" curvature %+.2f\n",
out[i].lambda, out[i].resid, out[i].curv);
}
std::printf("pool of %zu candidates: %d infeasible, %d not KKT, "
"%d KKT saddles,\n",
cand.size(), counts[0].load(), counts[1].load(),
counts[2].load());
std::printf(" %d strict local minimizers\n", counts[3].load());
}
Output:
(+1.00, +0.00): KKT + SOSC: strict local minimizer
lambda 0.500 residual 0.0e+00 curvature +1.00
(-1.00, +0.00): KKT, not a local minimizer
lambda 1.500 residual 0.0e+00 curvature -1.00
(+0.60, +0.80): feasible, not KKT
lambda 0.700 residual 8.0e-01 curvature +0.00
(+0.50, +0.00): infeasible
lambda 0.000 residual 0.0e+00 curvature +0.00
(+2.00, +0.00): feasible, not KKT
lambda 0.000 residual 3.0e+00 curvature +0.00
pool of 2005 candidates: 1 infeasible, 2002 not KKT, 1 KKT saddles,
1 strict local minimizers
The cost per candidate is a few dozen floating-point operations with no data dependence between candidates, and the only shared state is the four counters, which are atomic increments. That is the shape of work a GPU does well, and it is the shape of the incumbent side of a global solver: many points, one test each, verdicts that do not depend on the order of evaluation. On the circle the stationarity residual is \(|\sin t|\) and the least-squares multiplier is \(1 - \cos(t)/2\), so loosening \(\varepsilon_{\mathrm{stat}}\) lets a band of points around each KKT point pass as \(\varepsilon\)-KKT points, which is what a looser termination tolerance in a solver does.
The C++ check as a batch: 2,005 candidates, T threads, 4 counters
the pool 5 chosen by hand, 2,000 spread around the circle
| | | | | | | | |
v v v v v v v v v
threads thread t of T classifies the candidates t, t + T,
t + 2T, ...: feasibility, the multiplier, the
residual, the curvature; no data dependence
\__________________________ _______________________/
v
counters counts[code] += 1, an atomic increment: the only
shared state
|
v
infeasible not KKT KKT saddle strict local min
1 2002 1 1
Where this is used
Every global solver in Section 5 uses a local method as its supplier of incumbents and treats its output as a candidate. BARON runs a randomized multistart local search in its preprocessor, controlled by its NumLoc option, and launches further local searches at nodes of the tree. SCIP solves the NLP that remains after fixing the integer variables, dives with an NLP, that is, fixes variables one at a time and re-solves (Section 3.2), and uses the local solver inside the Undercover heuristic, a procedure for finding feasible points (Section 3.2). Couenne calls IPOPT for the upper bounds. Gurobi 13 added a nonlinear barrier method whose statuses are LOCALLY_OPTIMAL and LOCALLY_INFEASIBLE, and uses the same code inside its global solver to find feasible points.BARON, BARON User Manual and Installation Guide, v. 2026.9.10, The Optimization Firm, minlp.com/baron-user-manual, options NumLoc, DoLocal and NLPSol; S. Vigerske and A. Gleixner, "SCIP: global optimization of mixed-integer nonlinear programs in a branch-and-cut framework", Optimization Methods and Software 33 (2018); T. Berthold and A. M. Gleixner, "Undercover: a primal MINLP heuristic exploring a largest sub-MIP", Mathematical Programming 144 (2014); P. Belotti, J. Lee, L. Liberti, F. Margot and A. Wächter, "Branching and bounds tightening techniques for non-convex MINLP", Optimization Methods and Software 24 (2009); Gurobi Optimizer Reference Manual, version 13.0, parameters NLPHeur, NLBarPFeasTol, NLBarDFeasTol and NLBarCFeasTol, docs.gurobi.com. The two statuses, the preview-feature label and the caveat that the nonlinear barrier "does not guarantee a globally optimal solution, unless the problem is convex" are from the Gurobi 13.0 release notes, "Additions, changes and removals in Gurobi 13.0" (November 2025), docs.gurobi.com. In convex MINLP the local method is the engine rather than a heuristic. With the integer variables fixed, the remaining NLP is convex, its KKT point is its global optimum by Proposition 1.3.5, and the linearization at that point is a valid cut. That is the outer approximation method of Section 3.4. For a nonconvex MINLP the same subproblem is nonconvex, the local solver returns a KKT point that need not be the subproblem's optimum, and the linearization is no longer valid. That single difference is where the convex and the nonconvex methods of this series part ways.
What parallelizes
The parallel work divides along the same line as the certificates. Multistart and the verification of Algorithm 1.3.10 are independent across starts and candidates, and a GPU can run thousands of small local solves in lockstep on the same instructions and different data. Inside one local solve the function and derivative evaluations parallelize per constraint, and the dense fronts of a factorization parallelize within a front, but the pivot order of the indefinite factorization and its inertia test do not. That factorization is the reason a general nonlinear interior-point method does not run on a GPU as it stands, and Section 7.5 describes what has been done about it. The augmented Lagrangian family, which needs only products, reductions and projections, is the local method that fits the hardware without such surgery, and it is the ancestor of the ADMM codes (the alternating direction method of multipliers, Section 7.6) and the low-rank semidefinite codes of Section 7. The first-order linear programming codes of Section 7 come from a different lineage, the primal-dual methods of Section 7.2.
Integrality is a curve
Section 1.2 presented the two sources of nonconvexity in a MINLP as two kinds of hole: the set \(\{0, 1\}\) has a hole in the middle, and a curve that bends the wrong way has a hole beside it. This subsection shows that the first hole is a special case of the second. The constraint "\(x\) is \(0\) or \(1\)" is a concave quadratic inequality. Dropping it is the same operation as replacing that quadratic by its convex envelope, and the LP relaxation of an integer program is a chord relaxation, the term Proposition 3.5.11 fixes for replacing a concave term by the chord of its graph, which Section 2.4 shows is its convex envelope. The identification matters now for two reasons. It means that the relaxation machinery of Section 2 and the search of Section 3 treat integer and continuous nonconvexity by one mechanism. And it gives the smallest hard nonconvex program: a quadratic with one negative eigenvalue.
Two points as a parabola
The convex envelope of a function \(f\) on a box \(B\), written \(\operatorname{vex}_B f\), is the pointwise largest convex function that lies below \(f\) on \(B\). Section 2.4 studies it in general, and here only its value for one parabola is needed.
Proposition 1.4.1 (the integrality constraint as a concave quadratic). Let \(0 \le x \le 1\) and \(f(x) = x - x^2\). Then (i) \(x \in \{0, 1\}\) if and only if \(f(x) \le 0\); (ii) \(f(0) = f(1) = 0\), \(f > 0\) on \((0, 1)\), \(f''(x) = -2\), and \(\max_{[0,1]} f = \tfrac14\) at \(x = \tfrac12\); (iii) the convex envelope of \(f\) on \([0, 1]\) is the zero function; (iv) the relaxed constraint \(\operatorname{vex}_{[0,1]} f\,(x) \le 0\) has feasible set \([0, 1]\).
Proof. (i) and (ii): \(f(x) = x(1 - x)\) is the product of two nonnegative factors on \([0, 1]\), so \(f \ge 0\) there, with equality exactly when a factor vanishes. The derivative \(1 - 2x\) vanishes at \(\tfrac12\), where \(f = \tfrac14\). (iii): the zero function is convex and lies below \(f\) on \([0, 1]\), so \(\operatorname{vex}_{[0,1]} f \ge 0\). If \(h\) is convex with \(h \le f\) then \(h(0) \le 0\) and \(h(1) \le 0\), hence \(h(x) \le (1 - x) h(0) + x h(1) \le 0\) for every \(x \in [0, 1]\), so \(\operatorname{vex}_{[0,1]} f \le 0\). (iv) is immediate from (iii). ∎
(Two routes to one polytope) Read the four parts together. The integer program \(\min\{ c^\top x : Ax \le b,\ x \in \{0,1\}^n \}\) is the continuous program \(\min\{ c^\top x : Ax \le b,\ x_j - x_j^2 \le 0,\ 0 \le x_j \le 1 \}\) with \(n\) concave quadratic constraints. Relaxing each concave constraint by its convex envelope, which is what Section 2.4 does to every nonconvex term, replaces \(x_j - x_j^2 \le 0\) by \(0 \le 0\) and leaves the polytope \(\{Ax \le b,\ 0 \le x \le 1\}\). That polytope is the LP relaxation. Dropping integrality and relaxing the concave term are the same operation. The envelope sags \(\tfrac14\) under the curve at the midpoint, which Section 2.4 measures as the maximum error of an envelope, and that sag is the whole reason the relaxation of an integer program can return \(x_j = \tfrac12\).
Proposition 1.4.1 across n variables: two routes to one polytope
min{ c'x : Ax <= b, x in {0, 1}^n } the integer program
| |
| (i): x_j in {0, 1} if and only if |
| x_j - x_j^2 <= 0 on [0, 1] |
v |
min{ c'x : Ax <= b, x_j - x_j^2 <= 0, | drop the
0 <= x_j <= 1 } | integrality
| n concave quadratic constraints |
| |
| (iii): each x_j - x_j^2 replaced by its |
| convex envelope on [0, 1], the zero |
| function, that is by 0 <= 0 |
v v
min{ c'x : Ax <= b, 0 <= x <= 1 } the LP relaxation
one polytope at the end of both routes; the envelope sags 1/4
below x_j - x_j^2 at x_j = 1/2, which is why the relaxation can
return x_j = 1/2
(General integers) The same reading extends to general integers. For an integer variable with integer bounds \(l \le y \le u\), define \(g(t) = \operatorname{frac}(t)\,(1 - \operatorname{frac}(t))\) with \(\operatorname{frac}(t) = t - \lfloor t \rfloor\). On each unit cell \([k, k + 1]\) this is \((t - k)(k + 1 - t)\), the parabola of Proposition 1.4.1 translated, so \(g \ge 0\), \(g\) vanishes exactly at the integers, \(g'' = -2\) inside each cell and \(g = \tfrac14\) at each half-integer. The function has a convex kink at every integer, so it is piecewise concave rather than concave. The envelope argument of the proposition used only \(g \ge 0\) and \(g(l) = g(u) = 0\), and it gives \(\operatorname{vex}_{[l,u]} g = 0\) again. An integer variable \(y \in [l, u]\) is therefore the constraint \(g(y) \le 0\), and its relaxation is the interval.
Two other encodings of integrality appear in Section 1.5. The polynomial equation \(\prod_{k=l}^{u} (y - k) = 0\) has degree \(u - l + 1\), which grows with the range. The equation \(\sin(\pi y) = 0\) has no degree at all, and it is what carries the undecidability results of that subsection over to continuous problems with trigonometric functions.
The running example R1 of the plane figures (Section 1.6), the two-variable program over the polygon with six half-planes, is in this reading the continuous program with the constraints \(g(x) \le 0\) and \(g(y) \le 0\) added to the six linear ones. Its feasible set is the \(22\) lattice points inside the polygon. Replacing the two \(g\) constraints by \(0 \le 0\) gives back exactly the polygon, with its \(6\) vertices and an area of \(18.245\). The figure below draws both pictures. On the curve panel the probe reports \(x - x^2\) at the point you choose: \(0.21\) at the default \(x = 0.30\), where the constraint is violated and the relaxed constraint \(0 \le 0\) holds; the plane panel marks every point whose fractional part in \(x\) is the probe's, since \(g\) takes the same value on all of them. Its stats are \(f'' = -2\), one negative eigenvalue, the sag \(0.25\) at \(x = 0.5\), the two feasible points of the constraint, the \(22\) integer points of the polygon and its area, printed as \(18.25\). The plane shows the chain of parabolas \(g\) along each axis, and the thing to look for is their zero sets: the bumps vanish exactly on the lattice lines, which is where the integer points lie.
The following program checks the numbers the figure prints: the two feasible points, the second derivative and the sag, the polygon's vertices and area, and the count of lattice points when the integrality curve is applied to both coordinates.
# Integrality as a curve: the numbers fig-parabola prints.
#
# x - x^2 <= 0 on [0, 1], its chord envelope, and the running example's
# polygon as the same relaxation in two variables.
import numpy as np
f = lambda x: x - x**2
xs = np.linspace(0, 1, 100001)
feas = xs[f(xs) <= 1e-12]
print(f"feasible points of x - x^2 <= 0 on [0, 1]: "
f"{sorted(set(np.round(feas, 6)))} (2 points)")
print(f"f'' = {-2:+d} (one negative eigenvalue);")
print(f" max of x - x^2 = {f(xs).max():.4f} "
f"at x = {xs[f(xs).argmax()]:.2f}: the sag over the chord")
print("relaxed constraint 0 <= 0 on [0, 1]: feasible set [0, 1],")
print(f" an interval of length {xs[-1] - xs[0]:.0f}")
# the running example R1: six half-planes a.x + b.y <= r; the polygon is
# their intersection
H = [(2, 5, 24.5), (5, 2, 30.5), (-3, 4, 11), (1, -2, 4.2),
(-1, 0, 0), (0, -1, 0)]
def clip(poly, a, b, r):
"""Sutherland-Hodgman: clip a polygon against a x + b y <= r."""
out = []
for i in range(len(poly)):
P, Q = poly[i], poly[(i + 1) % len(poly)]
sP, sQ = a * P[0] + b * P[1] - r, a * Q[0] + b * Q[1] - r
if sP <= 1e-12:
out.append(P)
if ((sP < -1e-12 < sQ) or (sQ < -1e-12 < sP)
or (sP <= 1e-12 < sQ) or (sQ <= 1e-12 < sP)):
t = sP / (sP - sQ)
out.append((P[0] + t * (Q[0] - P[0]),
P[1] + t * (Q[1] - P[1])))
return out
poly = [(-1, -1), (9, -1), (9, 8), (-1, 8)]
for a, b, r in H:
poly = clip(poly, a, b, r)
area = 0.5 * abs(sum(poly[i][0] * poly[(i + 1) % len(poly)][1]
- poly[(i + 1) % len(poly)][0] * poly[i][1]
for i in range(len(poly))))
inside = lambda x, y: all(a * x + b * y <= r + 1e-9 for a, b, r in H)
# the integrality curve on every unit cell
g = lambda t: (t - np.floor(t)) * (1 - (t - np.floor(t)))
lattice = [(x, y) for x in range(0, 8) for y in range(0, 7)
if inside(x, y) and g(x) <= 1e-12 and g(y) <= 1e-12]
vertices = [(round(abs(p[0]), 3), round(abs(p[1]), 3)) for p in poly]
print(f"polygon: {len(poly)} vertices "
f"[{', '.join(str(v) for v in vertices[:3])},")
print(f" {', '.join(str(v) for v in vertices[3:])}], "
f"area {area:.3f}")
print(f"integer points in the polygon: {len(lattice)};")
print(f" g(t) = frac(t)(1 - frac(t)) peaks at {g(0.5):.2f} "
"on every cell")
Output:
feasible points of x - x^2 <= 0 on [0, 1]: [0.0, 1.0] (2 points)
f'' = -2 (one negative eigenvalue);
max of x - x^2 = 0.2500 at x = 0.50: the sag over the chord
relaxed constraint 0 <= 0 on [0, 1]: feasible set [0, 1],
an interval of length 1
polygon: 6 vertices [(4.2, 0.0), (5.783, 0.792), (4.929, 2.929),
(1.87, 4.152), (0.0, 2.75), (0.0, 0.0)], area 18.245
integer points in the polygon: 22;
g(t) = frac(t)(1 - frac(t)) peaks at 0.25 on every cell
The computation is a polygon clipping and a scan of \(56\) lattice points, and the scan is independent per point. Its point is that the integer feasible set and the relaxed feasible set are produced by the same code path with one constraint swapped for its envelope.
A concave quadratic encodes integrality
The identification runs the other way as well, and the converse direction is the one with consequences for hardness. If integrality is a concave quadratic, then a concave quadratic can encode integrality, and minimizing one over a polytope inherits the hardness of integer programming. The precise statement is sharper: one concave direction suffices.
Theorem 1.4.2 (nonconvex quadratic programming; Sahni 1974; Pardalos and Vavasis 1991; Vavasis 1990). (i) Minimizing a quadratic function \(x^\top H x + c^\top x\) over a polyhedron \(\{Ax \le b\}\) with rational data is NP-hard. (ii) It remains NP-hard when \(H\) has exactly one negative eigenvalue and all other eigenvalues equal to zero, that is, when the objective is concave along one direction and linear along every other. (iii) The decision version, "is there a feasible \(x\) with objective at most \(k\)", is in NP: whenever the answer is yes, a feasible point of polynomial bit size attains it.S. Sahni, "Computationally related problems", SIAM Journal on Computing 3 (1974); P. M. Pardalos and S. A. Vavasis, "Quadratic programming with one negative eigenvalue is NP-hard", Journal of Global Optimization 1 (1991), whose abstract states the result for "a concave quadratic function with one concave direction"; S. A. Vavasis, "Quadratic programming is in NP", Information Processing Letters 36 (1990). Part (ii) is proved in the paper; the reduction is not reproduced here.
Proof sketch of (i). Let \(P = \{x : Ax \le b,\ 0 \le x \le 1\}\) be a polytope whose integer points are to be decided, which is NP-complete by Karp's theorem (Theorem 1.5.2 below). By Proposition 1.4.1 the concave quadratic \(q(x) = \sum_j (x_j - x_j^2)\) is nonnegative on \(P\) and vanishes exactly at its \(0\)–\(1\) points, so \(\min_P q = 0\) if and only if \(P\) contains an integer point. The Hessian of \(q\) is \(-2I\), with \(n\) negative eigenvalues. Part (ii) says that the number of negative eigenvalues can be brought down to one without losing hardness, and that reduction is Pardalos and Vavasis's. Part (iii) is not obvious, because the minimizer of a quadratic may be irrational. Vavasis shows that a rational point of polynomial size with the required objective value exists nearby. ∎
Theorem 1.4.2 (i) and (ii): integer points via a concave q
does P = {x : Ax <= b, 0 <= x <= 1} contain an integer point?
NP-complete by Karp's theorem (Theorem 1.5.2)
|
| q(x) = sum_j (x_j - x_j^2): by Proposition 1.4.1,
| q >= 0 on P and q = 0 exactly at its 0-1 points
v
min_P q = 0 if and only if P contains an integer point;
the Hessian of q is -2I: n negative eigenvalues
|
| Pardalos and Vavasis: brought down to one without
| losing hardness (the reduction is not reproduced)
v
part (ii): one negative eigenvalue, all others zero, and
minimizing over a polyhedron is still NP-hard
(One concave direction suffices) Part (ii) says that hardness does not require many concave directions. A quadratic that curves the wrong way along a single direction, and is flat along all others, already forces a search over the vertices of a polytope, and the vertices of a polytope can encode a combinatorial problem. This is the lower limit of how nonconvex a problem must be to be hard, and it is the reason the two bilinear terms of Haverly's pooling problem, one per outflow of its single pool (R4 in Section 1.6, solved in Section 3.5), already make it a global optimization problem. For a fixed number \(t\) of negative eigenvalues, Vavasis later showed that an \(\varepsilon\)-approximate solution can be computed in time polynomial in \(n\) and \(1/\varepsilon\). A bounded number of concave directions therefore leaves the problem approximable, and Section 1.5 shows that without such a bound it is hard even to approximate.S. A. Vavasis, "Approximation algorithms for indefinite quadratic programming", Mathematical Programming 57 (1992).
The cleanest instance of a curve encoding a combinatorial quantity is older than any of this, and its proof fits in a paragraph. It concerns the clique number of a graph, the size of its largest clique, a set of vertices every two of which are joined by an edge, and computing that number is NP-hard.
Theorem 1.4.3 (Motzkin and Straus 1965). Let \(G\) be a graph with adjacency matrix \(A_G\) and clique number \(\omega(G)\), and let \(\Delta = \{ x \ge 0 : \sum_i x_i = 1 \}\) be the standard simplex. Then \(\max_{x \in \Delta} x^\top A_G\, x = 1 - 1/\omega(G)\).T. S. Motzkin and E. G. Straus, "Maxima for graphs and a new proof of a theorem of Turán", Canadian Journal of Mathematics 17 (1965).
Proof. For a clique \(K\) of size \(\omega\), the vector with weight \(1/\omega\) on each vertex of \(K\) and zero elsewhere gives \(x^\top A_G x = \omega(\omega - 1)/\omega^2 = 1 - 1/\omega\), so the maximum is at least \(1 - 1/\omega(G)\). Conversely take a maximizer \(x\). If two nonadjacent vertices \(i\) and \(j\) both carry positive weight, the objective contains no \(x_i x_j\) term. As a function of the transfer of weight from \(i\) to \(j\) along the simplex it is therefore linear, and moving all of the weight to one of the two ends does not decrease it. Repeating, there is a maximizer supported on a clique, of some size \(k \le \omega(G)\), and on a clique of size \(k\) the quadratic \(\sum_{i \ne j} x_i x_j = 1 - \sum_i x_i^2\) is at most \(1 - 1/k \le 1 - 1/\omega(G)\) by the Cauchy–Schwarz inequality. ∎
(A quadratic that counts cliques) The feasible set here is the simplest polytope there is, and the objective is a quadratic. All of the difficulty is in the sign of the quadratic: maximizing \(x^\top A_G x\) over the simplex computes the clique number, which is NP-hard to compute, so minimizing the nonconvex quadratic \(-x^\top A_G x\) over the simplex is NP-hard. Section 4.7 will meet this problem again as the standard quadratic program, whose exact convex reformulation lives in the copositive cone. For now it is the companion of Proposition 1.4.1. Proposition 1.4.1 writes an integrality constraint as a concave quadratic, and Theorem 1.4.3 shows that minimizing a nonconvex quadratic over a simplex computes a combinatorial quantity. Integer and continuous nonconvexity are therefore handled by one mechanism in Sections 2 and 3.
Where this is used
(What solvers take from the identification) No solver writes a binary variable as \(x - x^2 \le 0\). The practical traffic runs the other way. A solver that knows a variable is binary has a finite branching tree, cutting planes and propagation rules, which tighten one variable's bounds from the others' (Section 2.6), that a solver facing an anonymous concave constraint does not, so the integer structure is kept explicit wherever the model provides it. Section 4.7 shows how the lifted products of the reformulation-linearization technique use the identity \(x_j^2 = x_j\) of binaries as a linear fact. What the identification buys is understanding of the machinery. The LP relaxation of every MILP in Section 3 is the chord relaxation of this subsection. The integer branch \(y \le \lfloor v \rfloor\) or \(y \ge \lceil v \rceil\) of Section 3.1 is a spatial split of the interval of a concave term at a point where the term vanishes. The spatial branch of Section 3.5 is the same split applied where it does not.
What parallelizes
For the GPU programme the identification says that a bounding kernel built for convex envelopes of continuous terms serves both kinds of nonconvexity. The integer kind is the easier one, because its tree is finite and its envelope is exact at the points that matter.
What is known about hardness
The integers make a MINLP combinatorial: with \(p\) binary variables there are \(2^p\) patterns. Integer linear programming is NP-hard (Theorem 1.5.2 below), and so is anything that contains it. The curves make it nonconvex. The relaxation that drops integrality no longer solves to a certifiable optimum on its own, because the continuous problem can itself have many KKT points, as Example 1.3.6 showed, and a method must search those too. The two interact. A MILP's relaxation is an LP, solved exactly and, at the sizes in this series, in a fraction of a second. A nonconvex MINLP's relaxation is either an NLP whose answer is only local or a convex program that is only an approximation, and the search must branch on the continuous variables as well as the integer ones. This subsection collects what is actually known, and in what form:
- what NP-hardness says and does not say;
- which versions of the problem have no algorithm at all;
- why finite bounds on the variables matter, for three different reasons;
- what fixed dimension buys, and for which classes;
- why approximation does not rescue the nonconvex case;
- why "convex" is a declaration rather than a discovery.
A table at the end puts the results side by side with the four-letter table of Section 1.1. The subsection closes with the statement that makes solvers terminate at all: they certify \(\varepsilon\)-optimality, not optimality.
Vocabulary
In this subsection \(n\) is the total number of variables, \(p\) of them integer, in the single-vector notation of the Notation paragraph of the introduction.
Definition 1.5.1 (complexity vocabulary). The input size of an instance with rational data is the total number of bits of its coefficients. A decision problem is in P if some algorithm answers every instance in time polynomial in the input size. It is in NP if every "yes" instance has a certificate that can be checked in polynomial time, and in co-NP if every "no" instance has one. A problem is NP-hard if every problem in NP reduces to it in polynomial time, and NP-complete if it is also in NP. A problem is decidable if some algorithm answers every instance in finite time. Fixed dimension means that the number of variables \(n\) is a constant, so that a running time of the form \(n^{O(n)} \cdot \mathrm{poly}(\text{input size})\) counts as polynomial. The existential theory of the reals, \(\exists\mathbb{R}\), is the class of problems polynomially equivalent to deciding whether a system of polynomial equations and inequalities with integer coefficients has a real solution. It satisfies \(\mathrm{NP} \subseteq \exists\mathbb{R} \subseteq \mathrm{PSPACE}\), where PSPACE is the class of problems solvable in polynomial space.M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness (W. H. Freeman, 1979), for the framework; J. Canny, "Some algebraic and geometric computations in PSPACE", STOC 1988, for the PSPACE bound; M. Schaefer and D. Štefankovič, "Fixed points, Nash equilibria, and the existential theory of the reals", Theory of Computing Systems 60 (2017), for the class \(\exists\mathbb{R}\) as it is now used.
NP-hardness, and what it does not say
Theorem 1.5.2 (Karp 1972). Deciding whether \(\{ x \in \{0, 1\}^n : Ax \le b \}\) is nonempty, for integer \(A\) and \(b\), is NP-complete. Consequently MILP, convex MINLP and nonconvex MINLP are NP-hard, since each contains this problem as a special case.R. M. Karp, "Reducibility among combinatorial problems", in Complexity of Computer Computations (Plenum, 1972), where 0–1 integer programming is one of the twenty-one problems.
Proof sketch. Membership in NP: a feasible \(x\) is a certificate of \(n\) bits, and checking \(Ax \le b\) is polynomial. Hardness: a 3-SAT formula is a conjunction of clauses, each the disjunction of three literals, where a literal is a Boolean variable or its negation, and deciding whether such a formula is satisfiable is NP-complete. A clause \(\ell_1 \lor \ell_2 \lor \ell_3\) becomes the inequality \(\tilde\ell_1 + \tilde\ell_2 + \tilde\ell_3 \ge 1\), where \(\tilde\ell_k\) is \(x_j\) if the literal is the variable \(j\) and \(1 - x_j\) if it is its negation. The formula is satisfiable if and only if the system of clause inequalities has a \(0\)–\(1\) solution. ∎
Theorem 1.5.2: a 3-SAT formula as a system of 0-1 inequalities
formula clause and clause and ... and clause
| | |
v v v
system one inequality per clause, x in {0, 1}^n
clause l_1 or l_2 or l_3 ---> l~_1 + l~_2 + l~_3 >= 1
literal the variable j ---> l~_k = x_j
its negation ---> l~_k = 1 - x_j
satisfiable <==> the clause inequalities have a 0-1 solution
(What the theorem does not say) The lattice gives \(2^n\) candidates with no local structure to exploit, and linear inequalities are expressive enough to encode logic. Everything else in this subsection inherits at least this difficulty. Two things the theorem does not say should be said at once. It is a statement about the worst case over all instances of growing size: it says nothing about any particular instance, and nothing about instances drawn from a distribution. The simplex method is the standard example: Section 1.1 recorded that it is polynomial in the smoothed sense, and that is why it is fast in practice. Branch and bound solves random binary programs with a fixed number of constraints in a polynomial number of nodes with high probability. And yet there are explicit families, treated in Section 3.1, on which every branch-and-bound tree is exponential even after random perturbation of every coefficient.S. S. Dey, Y. Dubey and M. Molinaro, "Branch-and-bound solves random binary IPs in poly(n)-time", Mathematical Programming 200 (2023); S. S. Dey, Y. Dubey and M. Molinaro, "Lower bounds on the size of general branch-and-bound trees", Mathematical Programming 198 (2023), whose Theorem 2 is the perturbed family. The second thing is that hardness is not undecidability. The next theorem says that the decision version of integer linear programming, with unbounded integers, is in NP, so that an algorithm exists and a short certificate exists. Theorem 1.5.4 then says that quadratic constraints destroy both. The contrast with Theorem 1.5.7 below is worth holding in mind: a single quadratic inequality beside linear constraints keeps the problem in NP, while a system of quadratic equations in unbounded integers has no algorithm at all.
Theorem 1.5.3 (integer linear programming is in NP; Borosh and Treybig 1976; von zur Gathen and Sieveking 1978; Papadimitriou 1981). If \(\{ x \in \mathbb{Z}^n : Ax \le b \}\) with integer data is nonempty, it contains a point whose binary encoding has length polynomial in the input size. Hence the decision version of integer linear programming, with unbounded integer variables, is NP-complete.I. Borosh and L. B. Treybig, "Bounds on positive integral solutions of linear Diophantine equations", Proceedings of the American Mathematical Society 55 (1976); J. von zur Gathen and M. Sieveking, "A bound on solutions of linear integer equalities and inequalities", Proceedings of the American Mathematical Society 72 (1978); C. H. Papadimitriou, "On the complexity of integer programming", Journal of the ACM 28 (1981).
(Why unbounded integers are harmless for linear constraints) The rational polyhedron has vertices of polynomial size by Cramer's rule, and an integer point can be found within a bounded distance of a vertex. This is what makes unbounded integer variables harmless for linear constraints: the bounds that a solver infers for them are a convenience, not a necessity. The same statement fails as soon as the constraints are quadratic, and it fails in the strongest possible way.
Undecidability
Theorem 1.5.4 (Jeroslow 1973; after Matiyasevich 1970). There is no algorithm that, given integer data, solves the problem of minimizing a linear form over quadratic constraints in integer variables. The refinement of Jones (1982), as stated by Hemmecke, Köppe, Lee and Weismantel, is that there is no algorithm for minimizing a linear form over polynomial constraints in at most \(10\) integer variables.R. G. Jeroslow, "There cannot be any algorithm for integer programming with quadratic constraints", Operations Research 21 (1973). Yu. V. Matiyasevich, "Enumerable sets are Diophantine", Soviet Mathematics Doklady 11 (1970), for the negative solution of Hilbert's tenth problem; J. P. Jones, "Universal Diophantine equation", Journal of Symbolic Logic 47 (1982), for the bound on the number of unknowns; R. Hemmecke, M. Köppe, J. Lee and R. Weismantel, "Nonlinear integer programming", in 50 Years of Integer Programming 1958–2008 (Springer, 2010), Theorems 3 and 4, and M. Köppe, "On the complexity of nonlinear mixed-integer optimization", in Mixed Integer Nonlinear Programming, IMA Vol. 154 (Springer, 2012), Section 3, for the statements in the form used here.
Proof sketch. Hilbert's tenth problem asks for an algorithm that decides, for a polynomial \(P\) with integer coefficients, whether \(P(y) = 0\) has a solution in nonnegative integers. Matiyasevich, completing the work of Davis, Putnam and Robinson, proved that none exists. Now let \(P\) be given and consider the integer program
\[\min\ u \quad \text{subject to} \quad (1 - u)\, P(y) = 0, \qquad u \in \mathbb{Z}_{\ge 0},\ y \in \mathbb{Z}^k_{\ge 0} .\]If \(P\) has an integer zero \(y^\star\), then \((u, y) = (0, y^\star)\) is feasible and the optimal value is \(0\). If it has none, the constraint forces \(u = 1\), every \(y\) is then feasible, and the optimal value is \(1\). An algorithm for the integer program would therefore decide Hilbert's tenth problem. The constraint is a polynomial of high degree, but every product in it can be replaced by a new integer variable and a quadratic equation, \(w = y_i y_j\) or \(w = y_i w'\), until only quadratic equations remain. Each quadratic equation is a pair of quadratic inequalities. The linear objective and the quadratic constraints are Jeroslow's class. Jones's universal Diophantine equation needs only \(9\) unknowns, which with the variable \(u\) gives the bound of \(10\) for polynomial constraints. ∎
Example 1.5.5 (the gadget on two polynomials). Take \(P(x) = x^2 - 2\). Introducing \(w = x \cdot x\) and \(v = u \cdot w\), the program reads \(\min u\) subject to \(w - 2 - v + 2u = 0\), \(w = x^2\), \(v = u w\), with \(u, x, w, v\) nonnegative integers: two quadratic equations and one linear one. Since \(2\) is not a square, the optimal value is \(1\). Replace \(P\) by \(x^2 - 4\) and the optimal value is \(0\), attained at \(x = 2\), \(u = 0\). Solving every program of this shape would decide every Diophantine equation, so no algorithm solves every program of this shape.
Example 1.5.5: an equation P(x) = 0 becomes an integer program
does P(x) = x^2 - 2 = 0 have a solution in nonnegative integers?
|
| the program of Theorem 1.5.4
v
min u s.t. (1 - u) P(x) = 0, u, x nonnegative integers
|
| a new variable for each product:
v w = x * x and v = u * w
min u s.t. w - 2 - v + 2u = 0 one linear equation
w = x^2 two quadratic
v = u w equations
u, x, w, v nonnegative integers
|
v
2 is not a square, so u = 1 is forced: optimum 1
(replace P by x^2 - 4: optimum 0, attained at x = 2, u = 0)
(What the undecidability theorem covers) The series will go on to assume finite bounds everywhere, so what the theorem does and does not cover matters. Integer variables without bounds, together with the ability to multiply two variables, form a programming language for arithmetic, and arithmetic is undecidable. If every integer variable has finite bounds, the pure integer problem is decidable by enumeration, as Hemmecke and his co-authors remark: "for cases where finite bounds for all variables are known, an algorithm trivially exists". Decidable is not tractable, since the bounded problem is still NP-hard by Theorem 1.5.2.
(Continuous variables do not restore decidability) The mixed case adds a nuance that the pure case hides: continuous variables do not restore decidability. Del Pia, Dey and Molinaro record that feasibility of a system with \(3{,}424\) quadratic inequalities, \(58\) linear inequalities, \(58\) integer variables and \(1{,}711\) continuous variables is undecidable. The reduction is from a quartic equation in \(58\) nonnegative integer unknowns.A. Del Pia, S. S. Dey and M. Molinaro, "Mixed-integer quadratic programming is in NP", Mathematical Programming 162 (2017), Section 1.1.
(Bounded polynomial problems are decidable) With finite bounds on the integer variables, each integer assignment leaves a system of polynomial inequalities over the reals, and such systems are decidable. Tarski's quantifier elimination decides them, Renegar's algorithm does so in time \((sd)^{O(n)}\) for \(s\) polynomials of degree \(d\) in \(n\) variables, and Canny places the problem in PSPACE.J. Renegar, "On the computational complexity and geometry of the first-order theory of the reals, Part I", Journal of Symbolic Computation 13 (1992); J. Canny, STOC 1988, cited above; M. Schaefer and D. Štefankovič, Theory of Computing Systems 60 (2017), cited above, for the fact that the feasibility question for a nonconvex quadratically constrained program over the reals is complete for \(\exists\mathbb{R}\), so that whether it lies in NP is open. So a bounded MINLP with polynomial data is decidable, by enumeration of the integers and real algebra on each leaf.
(Transcendental functions) As soon as the functions include \(\exp\), \(\log\) or \(\sin\), even that is not known. Richardson proved that deciding whether an expression vanishes identically is undecidable for expressions built from rational constants, \(\pi\), one real variable, addition, multiplication, \(\sin\), \(\exp\) and the absolute value. Since \(y \in \mathbb{Z}\) if and only if \(\sin(\pi y) = 0\), a continuous program with trigonometric constraints and unbounded variables inherits Theorem 1.5.4 outright. A bounded integer variable \(l \le y \le u\) is instead encoded by the polynomial equation \(\prod_{k=l}^{u} (y - k) = 0\) of degree \(u - l + 1\), which is algebraic, and that is why bounded polynomial problems stay in the decidable case of the previous paragraph. Decidability of the real field with exponentiation alone is open and is known to follow from Schanuel's conjecture.D. Richardson, "Some undecidable problems involving elementary functions of a real variable", Journal of Symbolic Logic 33 (1968); A. J. Wilkie, "Schanuel's conjecture and the decidability of the real exponential field", in Algebraic Model Theory (Kluwer, 1997).
(What a solver actually does) None of this is what a solver does. No global solver runs quantifier elimination: they all compute an \(\varepsilon\)-optimal point on a compact box by branch and bound, which converges for continuous functions, and the exact decision problems of these paragraphs are never posed.R. Horst and H. Tuy, Global Optimization: Deterministic Approaches, 3rd ed. (Springer, 1996), Chapter IV, "Branch and Bound", pp. 115–178. The closing paragraphs of this subsection and Section 3.5 return to this.
Why bounds matter: three reasons
Every global solver in Section 5 needs finite bounds on the variables that appear in nonconvex terms, inferring them when it can, and uses them for far more than termination. Three distinct facts stand behind the requirement. Without bounds an infimum need not be attained. A feasible point can be too large to write down. And the second fact, in its strongest form, is what makes quadratic constraints with unbounded integers undecidable: no computable bound on the size of a solution exists.
Proposition 1.5.6 (existence of optima). (i) If the feasible set \(\mathcal F\) is nonempty and compact and \(f\) is continuous, the infimum \(z^\star\) is attained. In particular a MINLP with finite bounds on all variables and continuous \(f, g_i\) attains its infimum. (ii) For a MILP with rational data, \(\operatorname{conv}(\mathcal F)\) is a rational polyhedron, and a linear objective that is bounded below on \(\mathcal F \ne \emptyset\) attains its infimum, bounds or no bounds (Meyer 1974, stated and proved as Theorem 2.1.4). (iii) Without bounds on the integer variables a convex MINLP with rational data need not attain its infimum:
\[\inf\ \{ x : x\,y \ge 1,\ x \ge 0,\ y \in \mathbb{Z},\ y \ge 1 \} \;=\; 0 ,\]and no feasible point has \(x = 0\).R. R. Meyer, "On the existence of optimal solutions to integer and mixed-integer programming problems", Mathematical Programming 7 (1974). The example in (iii) is folklore and is attributed to no one here.
Proof. (i) \(\mathbb{R}^n \times \mathbb{Z}^p\) is closed, the box is compact and \(\{g \le 0\}\) is closed, so \(\mathcal F\) is compact, and a continuous function on a compact set attains its infimum. (ii) Proof in Meyer's paper. The argument is sketched at Theorem 2.1.4, which states the result in full. (iii) The feasible points are the pairs with \(x \ge 1/y\), so \(x = 1/y\) is feasible for every integer \(y \ge 1\) and the infimum is \(0\), while \(x = 0\) violates \(xy \ge 1\) for every \(y\). ∎
The constraint \(xy \ge 1\) on the positive quadrant is the convex side of the hyperbola of Proposition 1.2.5, a rotated second-order cone, so even the continuous relaxation of the example in (iii) is a convex program whose infimum is not attained. Add the bound \(y \le U\) and the optimum is \(1/U\) at \(y = U\): \(0.1\) for \(U = 10\), \(10^{-6}\) for \(U = 10^6\). The lattice is discrete, so with the continuous variables bounded the only way to lose the infimum is to run off to infinity along an integer coordinate while the objective keeps improving. A bound closes that road.
(The second reason: the size of a solution) The second reason is the size of the solution. Consider the ten inequalities \(x_1 \ge 2\) and \(x_{i+1} \ge x_i^2\) for \(i = 1, \dots, 9\), each convex and each a few bits long. The smallest feasible integer point is \(2, 4, 16, 256, 65{,}536, \dots\), with \(x_n = 2^{2^{n-1}}\), and \(x_{10} = 2^{512}\) has \(155\) decimal digits: every feasible point has a binary encoding exponential in the input size. The same phenomenon appears in two integer variables through Pell's equation \(x^2 - d\,y^2 = \pm 1\). The identity \(682^2 - 125 \cdot 61^2 = 465{,}124 - 465{,}125 = -1\) is a check the reader can do by hand that such equations have solutions at all. The fundamental solutions printed next show how fast the smallest solution grows with \(d\). For \(x^2 - 61 y^2 = 1\) it is \(x = 1{,}766{,}319{,}049\), ten digits for a two-digit \(d\), and for \(d = 109\) it has fifteen. Two quadratic inequalities in two integer variables can therefore force their smallest solution to be exponential in the bit size of the data. The following program prints the three computations.
# Why bounds matter: three computations.
#
# (a) An infimum that is not attained without a bound (Proposition
# 1.5.6 (iii)); (b) a feasible point whose size is doubly exponential;
# (c) Pell's equation, where the smallest solution is exponential in the
# input size.
from math import isqrt
print("(a) inf x s.t. x*y >= 1, x >= 0, y a positive integer:")
for U in (10, 1000, 10**6):
print(f" with y <= {U:>7}: optimum 1/U = {1/U:.6g} at y = {U};")
print(" without a bound the infimum 0 is never attained")
print("(b) x_1 >= 2, x_{i+1} >= x_i^2: the smallest feasible chain")
chain = [2]
for _ in range(9):
chain.append(chain[-1]**2)
print(f" first five: {chain[:5]};")
print(f" x_10 = 2^(2^9) has {len(str(chain[9]))} decimal digits")
print("(c) Pell: 682^2 - 125*61^2 =", 682**2 - 125*61**2)
def pell(d):
"""Fundamental solution of x^2 - d y^2 = 1, by continued fractions."""
a0 = isqrt(d)
m, q, a = 0, 1, a0
h1, h2, k1, k2 = 1, a0, 0, 1
while h2*h2 - d*k2*k2 != 1:
m = q*a - m
q = (d - m*m)//q
a = (a0 + m)//q
h1, h2, k1, k2 = h2, a*h2 + h1, k2, a*k2 + k1
return h2, k2
for d in (2, 13, 61, 109):
x, y = pell(d)
print(f" x^2 - {d:>3} y^2 = 1: smallest solution x = {x} "
f"({len(str(x))} digits),")
print(f" y = {y}; check {x*x - d*y*y}")
Output:
(a) inf x s.t. x*y >= 1, x >= 0, y a positive integer:
with y <= 10: optimum 1/U = 0.1 at y = 10;
without a bound the infimum 0 is never attained
with y <= 1000: optimum 1/U = 0.001 at y = 1000;
without a bound the infimum 0 is never attained
with y <= 1000000: optimum 1/U = 1e-06 at y = 1000000;
without a bound the infimum 0 is never attained
(b) x_1 >= 2, x_{i+1} >= x_i^2: the smallest feasible chain
first five: [2, 4, 16, 256, 65536];
x_10 = 2^(2^9) has 155 decimal digits
(c) Pell: 682^2 - 125*61^2 = -1
x^2 - 2 y^2 = 1: smallest solution x = 3 (1 digits),
y = 2; check 1
x^2 - 13 y^2 = 1: smallest solution x = 649 (3 digits),
y = 180; check 1
x^2 - 61 y^2 = 1: smallest solution x = 1766319049 (10 digits),
y = 226153980; check 1
x^2 - 109 y^2 = 1: smallest solution x = 158070671986249 (15 digits),
y = 15140424455100; check 1
The continued-fraction loop runs in a number of steps proportional to the length of the period of \(\sqrt d\), which is itself the quantity that grows. Nothing here parallelizes and nothing needs to. These two examples are exactly the ones Del Pia, Dey and Molinaro use to show that their theorem is tight.
The squaring chain x_1 >= 2, x_{i+1} >= x_i^2: its smallest point
x_1 x_2 x_3 x_4 x_5 x_10
2 -----> 4 -----> 16 ----> 256 ---> 65,536 ---> ... ---> 2^512
^2 ^2 ^2 ^2 ^2
each arrow squares: x_n = 2^(2^(n-1)), and x_10 has 155 decimal
digits, forced by ten inequalities each a few bits long
(The third reason: no computable bound at all) The third reason is the second one with every bound removed. Jeroslow's gadget in Theorem 1.5.4 needs exactly two ingredients: integer variables whose size no computable function bounds, and the product of two variables. The squaring chain shows the first ingredient at work with convex constraints alone. If a computable bound on the size of a smallest solution existed, the problem would be decidable by enumeration, so the undecidability of Theorem 1.5.4 is the statement that no such bound exists.
Theorem 1.5.7 (mixed-integer quadratic programming is in NP; Del Pia, Dey and Molinaro 2017). Let \(H\) be a rational symmetric matrix, \(c\), \(d\), \(A\) and \(b\) rational, and let \(\phi\) be the total bit size of the data. If the set
\[F \;=\; \{ x \in \mathbb{Z}^p \times \mathbb{R}^{n-p} : x^\top H x + c^\top x + d \le 0,\ Ax \le b \}\]is nonempty, it contains a point whose bit size is bounded by a polynomial in \(\phi\). Consequently the decision versions of integer and mixed-integer quadratic programming, with a nonconvex quadratic objective and linear constraints, are NP-complete.A. Del Pia, S. S. Dey and M. Molinaro, Mathematical Programming 162 (2017), Theorem 1 and Corollary 2; the squaring chain and the Pell example are in their Section 1.1. Proof in the paper.
(What the NP-membership theorem unifies) The single quadratic inequality is the objective-level cut "\(f(x) \le k\)", so this is the decision version of MIQP, and the theorem unifies Theorem 1.5.3 with Vavasis's result for continuous quadratic programs (Theorem 1.4.2(iii)). It is tight in the sense just computed: with two quadratic inequalities the smallest feasible integer point can have exponential size. Nothing like it is known for MIQCP, and whether feasibility of a mixed-integer quadratically constrained program is in NP is open. The authors also note that integer quadratic programming in two variables is polynomial and that its polynomiality in fixed dimension is open, which is the right moment to say what fixed dimension does buy.
Fixed dimension: Lenstra's theorem
For integer linear programming the parameter that governs the difficulty is the number of variables, not the size of the data. For every fixed \(n\) the problem is polynomial, and the proof is a geometric fact about convex bodies and lattices rather than a trick of encoding. A convex body that contains no lattice point is thin in some integer direction, and the thinness is bounded by a function of the dimension alone. Branching on that direction creates a number of children that does not depend on the data, and the recursion ends after \(n\) levels. Coordinate branching, which is what every production solver does, has no such guarantee, and the example at the end of this part shows the difference as a node count.
Definition 1.5.8 (lattices, width, flatness). A lattice in \(\mathbb{R}^n\) is a set \(L = \{ \sum_i m_i b_i : m \in \mathbb{Z}^n \}\) for linearly independent \(b_1, \dots, b_n\), its basis. Two bases of the same lattice differ by a unimodular matrix, an integer matrix of determinant \(\pm 1\), and the determinant \(\det L = |\det(b_1, \dots, b_n)|\) is independent of the basis. The standard lattice is \(\mathbb{Z}^n\). A vector \(c \in \mathbb{Z}^n \setminus \{0\}\) is primitive if the greatest common divisor of its entries is \(1\). For a convex body \(K\) (compact, convex, with nonempty interior) and \(c \ne 0\), the width of \(K\) in the direction \(c\) is
\[w_c(K) \;=\; \max_{x \in K} c^\top x \;-\; \min_{x \in K} c^\top x ,\]the length of the interval \(c^\top K\), and the lattice width is \(w(K) = \min\{ w_c(K) : c \in \mathbb{Z}^n \setminus \{0\} \}\). The minimum is attained at a primitive \(c\), since \(w_{dc} = d\, w_c\). \(K\) is lattice-point-free if \(K \cap \mathbb{Z}^n = \emptyset\). The flatness constant is \(\mathrm{Flat}(n) = \sup\{ w(K) : K \subset \mathbb{R}^n \text{ convex body with no lattice point in its interior} \}\). The supremum over lattice-point-free bodies is the same number, since shrinking a body slightly about an interior point removes its boundary lattice points and changes its width by as little as one likes.
The width is measured by the functional \(c^\top x\) and not by Euclidean distance. The reason is the following proposition: the lattice hyperplanes \(c^\top x = k\) are spaced one unit apart in the functional, so \(w_c(K)\) counts, up to rounding, how many of them meet \(K\).
Proposition 1.5.9 (lattice hyperplanes and slabs). Let \(c \in \mathbb{Z}^n\) be primitive. (i) \(\{ c^\top x : x \in \mathbb{Z}^n \} = \mathbb{Z}\), so every lattice point lies on one of the hyperplanes \(H_k = \{ c^\top x = k \}\), \(k \in \mathbb{Z}\), each of which contains lattice points, and consecutive hyperplanes are at Euclidean distance \(1/\|c\|_2\). (ii) A convex set \(K\) with \(c^\top K = [a, b]\) contains a lattice point only if \([a, b]\) contains an integer. In particular a lattice-point-free slab \(\{ a \le c^\top x \le b \}\) has \(b - a < 1\), and at most \(\lfloor w_c(K) \rfloor + 1\) of the hyperplanes \(H_k\) meet \(K\).
Proof. (i) \(c^\top x\) is an integer for integer \(x\), and by Bézout's identity there are integers \(m_i\) with \(\sum_i c_i m_i = \gcd(c) = 1\), so every integer \(k = c^\top (k m)\) is attained. (ii) If \(x \in K \cap \mathbb{Z}^n\) then \(c^\top x\) is an integer in \([a, b]\). An interval of length \(b - a \ge 1\) contains an integer, and an interval of length \(w\) contains at most \(\lfloor w \rfloor + 1\) of them. ∎
Theorem 1.5.10 (the flatness theorem; Khinchine 1948). For every \(n\) the constant \(\mathrm{Flat}(n)\) is finite: a convex body in \(\mathbb{R}^n\) whose interior contains no lattice point has lattice width at most \(\mathrm{Flat}(n)\).A. Khinchine, "A quantitative formulation of Kronecker's theory of approximation" (Russian), Izvestiya Akademii Nauk SSSR, Seriya Matematicheskaya 12 (1948), as cited for the theorem by W. Banaszczyk, A. E. Litvak, A. Pajor and S. J. Szarek, "The flatness theorem for nonsymmetric convex bodies via the local theory of Banach spaces", Mathematics of Operations Research 24 (1999). Known bounds: \(\mathrm{Flat}(n) = O(n^2)\) (R. Kannan and L. Lovász, "Covering minima and lattice-point-free convex bodies", Annals of Mathematics 128 (1988)); \(O(n^{3/2})\) (Banaszczyk, Litvak, Pajor and Szarek 1999); \(O(n^{4/3} \log^{O(1)} n)\) by combining that with M. Rudelson, "Distances between non-symmetric convex bodies and the \(MM^*\)-estimate", Positivity 4 (2000); \(O(n \log^{O(1)} n)\) (V. Reis and T. Rothvoss, "The subspace flatness conjecture and faster integer programming", FOCS 2023); \(\mathrm{Flat}(n) \ge n\) by the simplex \(\{x_i \ge \varepsilon, \sum_i x_i \le n - \varepsilon\}\), and \(\mathrm{Flat}(n) \ge 2n - o(n)\) (L. Mayrhofer, J. Schade and S. Weltge, "Lattice-free simplices with lattice width \(2d - o(d)\)", IPCO 2022, Lecture Notes in Computer Science 13265); \(\mathrm{Flat}(1) = 1\) and \(\mathrm{Flat}(2) = 1 + 2/\sqrt 3 \approx 2.15\) (C. A. J. Hurkens, "Blowing up convex sets in the plane", Linear Algebra and its Applications 134 (1990)), attained by a triangle with one lattice point on each edge; no other exact value is known, and \(\Theta(n)\) is conjectured.
Theorem 1.5.11 (integer programming in fixed dimension; Lenstra 1983). Fix \(n\). There is an algorithm that, given \(A \in \mathbb{Z}^{m \times n}\) and \(b \in \mathbb{Z}^m\) with \(m\) arbitrary, decides whether \(\{ x \in \mathbb{Z}^n : Ax \le b \}\) is nonempty, and finds a point of it if so, in time polynomial in the bit size of \((A, b)\). The same holds for minimizing a linear objective over that set, and for mixed-integer linear programs with \(n\) integer variables and any number of continuous ones.H. W. Lenstra, Jr., "Integer programming with a fixed number of variables", Mathematics of Operations Research 8 (1983). The dependence on \(n\) has been improved from Lenstra's \(2^{O(n^3)}\), or \(2^{O(n^2)}\) with the LLL-based count given below, to \(n^{O(n)}\) (R. Kannan, "Minkowski's convex body theorem and integer programming", Mathematics of Operations Research 12 (1987)), \(n^{4n/3 + o(n)}\) (D. Dadush, C. Peikert and S. Vempala, "Enumerative lattice algorithms in any norm via M-ellipsoid coverings", FOCS 2011), \(2^{O(n)} n^n\) (D. Dadush, PhD thesis, Georgia Institute of Technology, 2012) and \((\log 2n)^{O(n)}\) with randomization (Reis and Rothvoss, FOCS 2023); A. Frank and É. Tardos, "An application of simultaneous Diophantine approximation in combinatorial optimization", Combinatorica 7 (1987), make the running time strongly polynomial.
The mechanism of the proof is worth having in full. It is also the mechanism behind the convex case (Theorem 1.5.17), and it says precisely which part of the work depends on the data and which on the dimension. It needs two tools: an ellipsoidal rounding of the body and a reduced basis of the lattice.
Definition 1.5.12 (rounding; reduced bases). A \(\rho\)-rounding of a convex body \(K\) is an ellipsoid \(E = p + A B_2^n\), with \(A\) invertible and \(B_2^n\) the unit ball, such that \(E \subseteq K \subseteq p + \rho A B_2^n\). Every convex body has an \(n\)-rounding, its maximal-volume inscribed ellipsoid dilated by \(n\) (John's theorem), and Lenstra constructs one for a polytope given by inequalities in polynomial time with a \(\rho\) depending on \(n\) alone. A parallelotope \(p + A[-1, 1]^n\) has the \(\sqrt n\)-rounding \(p + A B_2^n\), since \(B_2^n \subseteq [-1, 1]^n \subseteq \sqrt n\, B_2^n\). For a lattice basis \(b_1, \dots, b_n\) the Gram–Schmidt vectors are \(b_1^* = b_1\) and \(b_i^* = b_i - \sum_{j < i} \mu_{ij} b_j^*\) with \(\mu_{ij} = b_i^\top b_j^* / \|b_j^*\|^2\), so that \(\det L = \prod_i \|b_i^*\|\). The basis is LLL-reduced (with parameter \(3/4\)) if \(|\mu_{ij}| \le \tfrac12\) for all \(j < i\) and \(\|b_i^*\|^2 \ge (\tfrac34 - \mu_{i,i-1}^2) \|b_{i-1}^*\|^2\) for \(i = 2, \dots, n\).A. K. Lenstra, H. W. Lenstra, Jr. and L. Lovász, "Factoring polynomials with rational coefficients", Mathematische Annalen 261 (1982), which gives the reduction algorithm and proves that it runs in time polynomial in \(n\) and in the bit size of the input basis.
Proposition 1.5.13 (what reduction buys). For an LLL-reduced basis: (i) \(\|b_i^*\|^2 \ge \tfrac12 \|b_{i-1}^*\|^2\), hence \(\|b_j^*\|^2 \le 2^{\,i-j} \|b_i^*\|^2\) for \(j \le i\); (ii) \(\|b_i\| \le 2^{(i-1)/2} \|b_i^*\|\); (iii) \(\max_i \|b_i\| \le 2^{(n-1)/2} \|b_n^*\|\); (iv) the lattice is contained in the union of the parallel hyperplanes \(k\, b_n + \operatorname{span}(b_1, \dots, b_{n-1})\), \(k \in \mathbb{Z}\), and consecutive ones are at Euclidean distance \(\|b_n^*\|\).
Proof. (i) is the Lovász condition with \(\mu_{i,i-1}^2 \le \tfrac14\), iterated. (ii): \(\|b_i\|^2 = \|b_i^*\|^2 + \sum_{j<i} \mu_{ij}^2 \|b_j^*\|^2 \le \|b_i^*\|^2 (1 + \tfrac14 \sum_{j<i} 2^{i-j}) \le 2^{i-1} \|b_i^*\|^2\). (iii) combines (ii) with (i) for \(j = i\), \(i = n\). (iv): a lattice vector \(\sum_i m_i b_i\) lies on the hyperplane with \(k = m_n\), and the component of \(b_n\) orthogonal to the span of the others is \(b_n^*\). ∎
Theorem 1.5.14 (flatness through rounding and reduction; the step of Lenstra's algorithm). Let \(K \subset \mathbb{R}^n\) be a convex body with a \(\rho\)-rounding \(E = p + A B_2^n\), let \(L = A^{-1} \mathbb{Z}^n\) with an LLL-reduced basis \(b_1, \dots, b_n\), write \(c = A^{-1} p = \sum_i \alpha_i b_i\) and let \(y = \sum_i \lfloor \alpha_i \rceil\, b_i\), rounding each coefficient to the nearest integer. Then at least one of the following holds: (a) \(A y \in K \cap \mathbb{Z}^n\); (b) there is a primitive \(w \in \mathbb{Z}^n\) with \(w_w(K) < \rho\, n\, 2^{(n-1)/2}\), so that at most \(\rho\, n\, 2^{(n-1)/2} + 1\) of the lattice hyperplanes \(w^\top x = k\) meet \(K\). In particular, if \(K\) is lattice-point-free then \(w(K) < \rho n 2^{(n-1)/2}\), and with John's \(\rho = n\) this gives \(\mathrm{Flat}(n) \le n^2 2^{(n-1)/2}\).
Proof. Let \(\tau(x) = A^{-1} x\), so that \(B(c, 1) \subseteq \tau(K) \subseteq B(c, \rho)\) and \(\tau(\mathbb{Z}^n) = L\). Since each \(|\alpha_i - \lfloor \alpha_i \rceil| \le \tfrac12\),
\[\|c - y\| \;=\; \Big\| \sum_i (\alpha_i - \lfloor \alpha_i \rceil)\, b_i \Big\| \;\le\; \tfrac12 \sum_i \|b_i\| \;\le\; \tfrac n2 \max_i \|b_i\| .\]If \(\max_i \|b_i\| \le 2/n\) then \(y \in B(c, 1) \subseteq \tau(K)\), so \(Ay \in K\), and \(Ay \in \mathbb{Z}^n\) because \(y \in L\): this is (a). Otherwise \(\max_i \|b_i\| > 2/n\) and Proposition 1.5.13(iii) gives \(\|b_n^*\| > 2^{-(n-1)/2} \cdot 2/n\). Let \(\varphi\) be the linear functional with \(\varphi(b_i) = 0\) for \(i < n\) and \(\varphi(b_n) = 1\), whose level sets are the hyperplanes of Proposition 1.5.13(iv), spaced \(\|b_n^*\|\) apart. It maps the ball \(B(c, \rho)\) onto an interval of length \(2\rho / \|b_n^*\| < \rho n 2^{(n-1)/2}\), and \(\varphi(\tau(K))\) is a subinterval. The functional \(\psi(x) = \varphi(A^{-1} x)\) takes integer values on \(\mathbb{Z}^n\), because \(\psi(\mathbb{Z}^n) = \varphi(L) = \mathbb{Z}\), so \(\psi(x) = w^\top x\) for the integer vector \(w = (\psi(e_1), \dots, \psi(e_n))\), and \(w\) is primitive because its values on the lattice are all of \(\mathbb{Z}\) rather than a proper subgroup. The width of \(K\) in the direction \(w\) is the length of \(\psi(K) = \varphi(\tau(K))\), which is less than \(\rho n 2^{(n-1)/2}\), and Proposition 1.5.9(ii) bounds the number of hyperplanes: this is (b). If \(K\) is lattice-point-free, (a) is impossible, so (b) holds and \(w(K) \le w_w(K)\). ∎
(The step in a picture, and its constants) In a picture: after the rounding the body is nearly a ball, and a reduced basis of the lattice is nearly orthogonal. Either every basis vector is short compared with the inscribed ball, in which case rounding the coordinates of the centre lands on a lattice point inside the body. Or some basis vector is long, in which case the last Gram–Schmidt vector is long, the lattice is a stack of widely spaced parallel hyperplanes, and a ball of radius \(\rho\) meets few of them. The constants are poor. For \(n = 2\) the theorem gives \(4\sqrt 2 \approx 5.66\) against the true \(\mathrm{Flat}(2) = 1 + 2/\sqrt 3 \approx 2.15\), and for \(n = 10\) it gives about \(2{,}263\). What matters is that the bound involves \(\rho\) and \(n\) and nothing from \(A\) or \(b\).
Theorem 1.5.14: the step of Lenstra's algorithm
K, with a rho-rounding E = p + A B_2^n
|
| tau(x) = A^{-1} x
v
B(c, 1) inside tau(K) inside B(c, rho), with c = A^{-1} p;
the lattice Z^n becomes L = A^{-1} Z^n
|
| an LLL-reduced basis b_1, ..., b_n of L
v
c = sum_i alpha_i b_i --> y = sum_i round(alpha_i) b_i,
||c - y|| <= (n/2) max_i ||b_i||
|
+-- max_i ||b_i|| <= 2/n --> (a) y in B(c, 1), so A y is
| a lattice point of K
|
+-- max_i ||b_i|| > 2/n ---> ||b_n*|| > 2^{-(n-1)/2} 2/n:
L lies on hyperplanes spaced
||b_n*|| apart, few of which
meet B(c, rho)
(b) a primitive w with
w_w(K) < rho n 2^{(n-1)/2}
Two terms in the algorithm box need a gloss. A separation oracle for a convex body is a routine that, given a point, either confirms that the point lies in the body or returns a hyperplane separating the point from it. The Hermite normal form of an integer matrix is the integer analogue of row echelon form, computable in polynomial time, from which the integer solutions of a system of linear equations are read off.
Algorithm 1.5.15 (Lenstra's algorithm for integer feasibility in fixed dimension).
Input a polytope P = {x in R^n : A x <= b} with integer data (or a
convex body given by a separation oracle and a bounding box).
Output a point of P ∩ Z^n, or the certificate "P ∩ Z^n is empty".
1. (empty) Solve an LP. If P is empty, return "empty".
2. (dimension) If P is not full-dimensional, find an equation
a^T x = beta valid on P (an LP). Its integer solutions
are empty or an affine lattice of dimension n - 1
(Hermite normal form); rewrite P in integer
coordinates on that lattice and recurse with n - 1.
3. (bounds) If P is unbounded, intersect it with a box of
polynomially bounded size that contains an integer
point of P whenever P has one (the size bound of
Theorem 1.5.3).
4. (rounding) Compute p and an invertible A with
p + A B ⊆ P ⊆ p + rho(n) A B,
rho(n) depending on n alone.
5. (reduction) LLL-reduce the basis A^{-1} e_1, ..., A^{-1} e_n of
L = A^{-1} Z^n; call the result b_1, ..., b_n.
6. (round) Write A^{-1} p = sum_i alpha_i b_i and
y = sum_i round(alpha_i) b_i.
If A y is in P, return A y.
7. (direction) Let psi(x) = coefficient of b_n in A^{-1} x = w^T x
with w in Z^n primitive (Theorem 1.5.14).
Compute k_min = ceil(min_P w^T x) and
k_max = floor(max_P w^T x) by two LPs.
[k_max - k_min + 1 <= rho(n) n 2^{(n-1)/2} + 1.]
8. (branch) For k = k_min, ..., k_max:
complete w to a unimodular U with first row w;
in the coordinates z = U x the hyperplane
w^T x = k is z_1 = k, so P ∩ {w^T x = k} is an
instance in the n - 1 integer variables
z_2, ..., z_n.
Recurse. If a call returns a point, map it back
by x = U^{-1} z and return it.
9. Return "empty".
Invariant
at every node P ∩ Z^n is the union over the children of
(child ∩ Z^n), so no integer point is lost; each child has
dimension n - 1; the fan-out of a node of dimension k is at most
F(k) = rho(k) k 2^{(k-1)/2} + 1;
the depth is at most n. A one-dimensional instance is an interval
and is solved by rounding its ends.
The tree of Algorithm 1.5.15: each level drops the dimension by one
dim n P
|
+----------------+---------------------+
| | |
dim n - 1 w^T x = k_min w^T x = k_min + 1 ... w^T x = k_max
/ | \
dim n - 2 . . .
:
dim 1 intervals, solved by rounding their ends
a child: P cap {w^T x = k} in the integer variables z_2, ..., z_n
(step 8), or P on an affine lattice of dimension n - 1 (step 2)
fan-out at dimension k at most F(k) = rho(k) k 2^{(k-1)/2} + 1,
no quantity from A or b; depth at most n; at most
n prod_{k=2}^{n} F(k) nodes, a constant for fixed n
(The cost of a node, and what parallelizes) The cost of a node is a constant number of linear programs over the current polytope, one rounding, one LLL reduction of an \(n \times n\) basis and one unimodular change of coordinates, all polynomial in the input size. The tree has at most \(n \prod_{k=2}^{n} F(k)\) nodes, a constant for fixed \(n\), which is \(2^{O(n^2)}\) for polynomial \(\rho\), and minimizing a linear objective is a binary search on its value with \(O(\text{input size})\) runs of the feasibility algorithm. What parallelizes is the children of a node, which are independent instances of dimension \(n - 1\) and can be solved concurrently with the first success cancelling its siblings, together with the linear programs of steps 1 to 3 and 7, which are the data-dependent part and can use any parallel LP method. The reduction of step 5 is a sequence of dependent swaps and stays sequential. The count depends on the dimension and not on the data because the fan-out bound \(F(k)\) contains \(\rho(k)\), \(k\) and \(2^{(k-1)/2}\) and no quantity from \(A\) or \(b\). The data enter only through the polynomial cost of each node and the sizes of the numbers carried along. Coordinate branching has no such bound, and the next example is the smallest one in which the difference is visible.
The strip: coordinate branching against the flat direction
The worked example of this part is the strip
\[K_N \;=\; \{ (x, y) \in \mathbb{R}^2 : 0 \le x \le N,\ 1 \le 10y - 5x \le 4 \}, \qquad N \ge 1 \text{ an integer},\]a parallelogram with vertices \((0, 0.1)\), \((0, 0.4)\), \((N, N/2 + 0.1)\) and \((N, N/2 + 0.4)\). Since \(10y - 5x = 5(2y - x)\), the two slanted inequalities say \(0.2 \le 2y - x \le 0.8\): the strip is a piece of a slab around the line \(2y - x = \tfrac12\), which is not a coordinate line. Consider the following search along a fixed integer form \(c^\top p\), which is Lenstra's step with the bookkeeping of a branch-and-bound tree. Branch and bound, defined in general in Section 3.1, splits the set to be searched into pieces and discards a piece as soon as it can be shown to contain no solution. Here each piece, called a node, is the strip cut down to a range of one integer form, and a piece is discarded, or closed (pruned, in the vocabulary of Section 3.1), when its relaxation is empty or when some coordinate's range holds no integer. In detail, a node carries integer bounds \(k_{\mathrm{lo}} \le c^\top p \le k_{\mathrm{hi}}\) and the relaxation \(K_N \cap \{ k_{\mathrm{lo}} \le c^\top p \le k_{\mathrm{hi}} \}\). The node closes if the relaxation is empty, or if the range of \(x\) or the range of \(y\) on the relaxation contains no integer, which is the bound rounding every solver applies (Section 2.6). Otherwise the range \([lo, hi]\) of \(c^\top p\) on the relaxation is computed. If it contains the integers \(k_1 < \dots < k_r\) the node gets the \(r\) children \(c^\top p = k_j\), and if it contains none it gets the two children \(c^\top p \le \lfloor lo \rfloor\) and \(c^\top p \ge \lceil hi \rceil\).
The search along an integer form c.p: one node
node: integer bounds k_lo <= c.p <= k_hi
relaxation: K_N cap {k_lo <= c.p <= k_hi}
|
relaxation empty? ---------------------- yes --> closed
| no
range of x or of y holds no integer? --- yes --> closed
| no (bound rounding)
[lo, hi] = the range of c.p on the relaxation
|
integers k_1 < ... < k_r in [lo, hi]? -- yes --> r children,
| no c.p = k_j
v
two children: c.p <= floor(lo) and c.p >= ceil(hi)
Proposition 1.5.16 (the strip). (i) \(K_N\) contains no lattice point: at \(x = k\) the admissible \(y\) form the interval \([k/2 + 0.1,\ k/2 + 0.4]\) of length \(0.3\), which is \([m + 0.1, m + 0.4]\) for \(k = 2m\) and \([m + 0.6, m + 0.9]\) for \(k = 2m + 1\). (ii) For a primitive direction \((a, b)\) the width is \(w_{(a,b)}(K_N) = |a + b/2|\, N + 0.3\,|b|\) when \(b \ne 0\) and \(|a| N\) when \(b = 0\). Hence \(w(K_N) = 0.6\), attained only along \(\pm(-1, 2)\), and every other primitive direction has width at least \(N/2 + 0.3\). (iii) For every integer \(N \ge 2\) the search above processes
\[N + 2 \ \text{ nodes along } x, \qquad \lfloor N/2 + 0.4 \rfloor + 1 \ \text{ along } y, \qquad \lfloor 1.5N + 0.4 \rfloor + 1 \ \text{ along } x + y, \qquad 3 \ \text{ along } 2y - x ,\]and no run finds a lattice point. For \(N = 1\) the root's range of \(y\) is \([0.1, 0.9]\), so every form closes the root at once.
Proof. (i) is read off \(y \in [x/2 + 0.1, x/2 + 0.4]\). Alternatively, \((-1, 2)\) is primitive and \(2y - x\) ranges over \([0.2, 0.8]\) on \(K_N\), so Proposition 1.5.9(ii) applies. (ii) The values of \(ax + by\) at the four vertices are \(0.1b\), \(0.4b\), \(s + 0.1b\) and \(s + 0.4b\) with \(s = (a + b/2) N\), and a linear functional on a polygon attains its extremes at vertices, so the width is \(|s| + 0.3|b|\). For \(b\) odd, \(|a + b/2| \ge \tfrac12\) and the width is at least \(N/2 + 0.3\), with equality for \(\pm(0, 1)\) and \(\pm(-1, 1)\). For \(b\) even and nonzero, either \(a + b/2 = 0\), which for a primitive vector means \((a, b) = \pm(-1, 2)\) and gives \(0.6\), or \(|a + b/2| \ge 1\) and the width is at least \(N + 0.6\). For \(b = 0\) the width is \(|a| N \ge N\). (iii) On \(K_N\) the ranges are \(x \in [0, N]\), \(y \in [0.1, N/2 + 0.4]\), \(x + y \in [0.1, 1.5N + 0.4]\) and \(2y - x \in [0.2, 0.8]\). For \(N \ge 2\) the ranges of \(x\) and of \(y\) contain integers, so the root does not close and branches on its form. Along \(x\) the children are \(x = k\) for \(k = 0, \dots, N\), and on each the range of \(y\) holds no integer by (i), so each closes: \(1 + (N + 1)\) nodes. Along \(y\) the children are \(y = k\) for \(k = 1, \dots, \lfloor N/2 + 0.4 \rfloor\), and on \(y = k\) the strip gives \(x \in [2k - 0.8, 2k - 0.2]\), no integer. Along \(x + y\) the children are \(x + y = k\) for \(k = 1, \dots, \lfloor 1.5N + 0.4 \rfloor\). Substituting \(x = k - y\) into \(0.2 \le 2y - x \le 0.8\) gives \(x \in [(2k - 0.8)/3, (2k - 0.2)/3]\), and an integer \(m\) in that interval would need \(3m - 2k \in [-0.8, -0.2]\), which no integer satisfies. Along \(2y - x\) the root range \([0.2, 0.8]\) holds no integer, so the root gets the children \(2y - x \le 0\) and \(2y - x \ge 1\), both with empty relaxation: \(3\) nodes. No child contains a lattice point because \(K_N\) contains none. ∎
The four searches on K_10 (Proposition 1.5.16 (iii), N = 10)
along x, 12 nodes: the root and the children x = k
root
|
+----+----+----+----+----+----+----+----+----+----+
| | | | | | | | | | |
0 1 2 3 4 5 6 7 8 9 10
along y, 6 nodes: the root and the children y = k
root
|
+----+----+----+----+
| | | | |
1 2 3 4 5
along x + y, 16 nodes: the root and the children x + y = k
root
|
+---+---+---+---+---+---+---+---+---+---+---+---+---+---+
| | | | | | | | | | | | | | |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
along 2y - x, 3 nodes: the root and one disjunction
root
|
+-----------+-----------+
| |
2y - x <= 0 2y - x >= 1
(empty) (empty)
every child is a leaf: along x its range of y holds no integer,
along y and x + y its range of x holds none, and along 2y - x
its relaxation is empty
The counts are the figure's. In its default view, \(N = 10\) with branching on \(x\), the search processes \(12\) nodes, the root and the eleven children \(x = 0, \dots, 10\), each closed because its range of \(y\) holds no integer, against \(3\) nodes along \(2y - x\). Along \(y\) it takes \(6\) and along \(x + y\) it takes \(16\). At the slider's end, \(N = 40\), the counts are \(42\), \(21\), \(61\) and \(3\). A hand count of the same search records only that the number of nodes along \(x\) is linear in \(N\), which is the same fact without the constant. The figure prints the exact \(N + 2\). Its lattice-width search over the primitive directions with \(|a|, |b| \le 6\) finds \(0.6\) along \((-1, 2)\) and then \(5.3\) along \((0, 1)\) and \((-1, 1)\) at \(N = 10\), which is part (ii) with \(N/2 + 0.3 = 5.3\). The Euclidean thickness of the strip in its thin direction is \(0.6/\sqrt 5 \approx 0.27\), while the lattice lines \(2y - x = k\) are \(1/\sqrt 5\) apart, so the strip fits between two consecutive lattice lines with room to spare. That is the picture to look for when the figure is set to \(2y - x\): the two half-planes hatched red and the strip in the white gap between them.
(Lenstra's step on the strip, by hand) Lenstra's step finds the thin direction by arithmetic, and on the strip with \(N = 10\) the arithmetic is short. The rounding of Definition 1.5.12 is the ellipse \(p + A B_2^2\) with \(p = (5, 2.75)\) and \(A = \big(\begin{smallmatrix} 5 & 0 \\ 2.5 & 0.15 \end{smallmatrix}\big)\), a \(\sqrt 2\)-rounding since \(K_{10} = p + A[-1, 1]^2\). Gauss reduction of the basis of \(L = A^{-1} \mathbb{Z}^2\) returns the short vector \(v_1 = (\tfrac25, 0)\), whose pull-back \(A v_1 = (2, 1)\) points along the level lines of \(2y - x\), and the lattice lines parallel to \(v_1\) are \(\tfrac{10}{3}\) apart in the transformed plane while the body lies in a disc of radius \(\sqrt 2\). No lattice line meets the body, so step 8 creates no children, step 9 returns "empty", and \(K_{10} \cap \mathbb{Z}^2 = \emptyset\) is proved in one step, along a direction of width \(0.6\) where Theorem 1.5.14 promises only a width below \(4\).The arithmetic in full: \(A^{-1} e_1 = (\tfrac15, -\tfrac{10}{3})\) and \(A^{-1} e_2 = (0, \tfrac{20}{3})\) span \(L\) with \(\det L = \tfrac43\). Gauss reduction, which is LLL in two dimensions, gives \(v_1 = 2 A^{-1} e_1 + A^{-1} e_2 = (\tfrac25, 0)\) and \(v_2 = A^{-1} e_1\), with \(\|v_1\| = 0.4\) and \(\|v_2^*\| = \det L / \|v_1\| = \tfrac{10}{3}\). The centre is \(c = A^{-1} p = (1, \tfrac53)\). The functional "coefficient of \(v_2\)" is \(\pm(2y - x)\), and the lattice points sit at second coordinate \(\tfrac{10}{3}(2y - x) - \tfrac53 \in \{ \pm\tfrac53, \pm 5, \dots \}\) relative to \(c\), while the body occupies \([-1, 1]\). The rounding test of step 6 fails first, as it must, because the nearest lattice point to \(c\) in the transformed metric, \(\tau(5, 3)\), is at distance \(\tfrac53 > 1\). The bound of Theorem 1.5.14 with \(n = 2\) and \(\rho = \sqrt 2\) is \(\rho n 2^{(n-1)/2} = 4\).
(The same arithmetic as preprocessing) The same arithmetic is the two-variable case of a preprocessing step. The matrix \(U = \big(\begin{smallmatrix} -1 & 2 \\ 0 & 1 \end{smallmatrix}\big)\) has determinant \(-1\), so \((s, t) = (2y - x, y)\) is a unimodular change of integer variables, with inverse \((x, y) = (2t - s, t)\). In the new variables the strip is \(\{ 0.2 \le s \le 0.8,\ 0 \le 2t - s \le N \}\). The integer variable \(s\) has the bounds \([0.2, 0.8]\), rounding them gives \(1 \le s \le 0\), and the root is infeasible by bound propagation (Section 2.6) alone, for every \(N\), with no branching. The lattice reformulation of Aardal, Hurkens and Lenstra for systems of linear Diophantine equations with bounds is this change of variables in general, computed by basis reduction, and it solved market-split instances that coordinate branching could not.K. Aardal, C. A. J. Hurkens and A. K. Lenstra, "Solving a system of linear Diophantine equations with lower and upper bounds on the variables", Mathematics of Operations Research 25 (2000), compute by lattice basis reduction a short solution \(x_0\) of \(Ax = d\) and a reduced basis \(B\) of the lattice \(\{ x \in \mathbb{Z}^n : Ax = 0 \}\), and branch on the coordinates \(\lambda\) of \(x = x_0 + B\lambda\). Because \(B\) is reduced, the reformulated polytope tends to be thin in its last coordinates, where ordinary bound rounding and variable branching then work. The market-split instances are those of G. Cornuéjols and M. Dawande, "A class of hard small 0-1 programs", INFORMS Journal on Computing 11 (1999), offered as a challenge because conventional branch and bound found even small ones extremely difficult. The reformulation solved instances with up to \(7\) equations and \(60\) variables, where the largest solved before had \(6\) equations and \(50\) variables for feasibility and \(4\) equations and \(30\) variables for optimization: K. Aardal, R. E. Bixby, C. A. J. Hurkens, A. K. Lenstra and J. W. Smeltink, "Market split and basis reduction: towards a solution of the Cornuéjols–Dawande instances", INFORMS Journal on Computing 12 (2000), abstract. K. Aardal and A. K. Lenstra, "Hard equality constrained integer knapsacks", Mathematics of Operations Research 29 (2004), explain the effect through the small number of lattice hyperplanes that meet the reformulated polytope in the direction of its last coordinate.
The unimodular change of variables on the strip
(s, t) = (2y - x, y) = U (x, y)
U = [ -1 2 ] det U = -1
[ 0 1 ]
(x, y) ------------------------------> (s, t)
<------------------------------
(x, y) = (2t - s, t)
0 <= x <= N <==> 0 <= 2t - s <= N
1 <= 10y - 5x <= 4 <==> 0.2 <= s <= 0.8
s is an integer variable: its bounds round to 1 <= s <= 0, so
the root is infeasible by bound propagation alone, for every N,
with no branching
The program below re-implements the figure's search and prints the node counts of Proposition 1.5.16 for the four forms. It checks the closed forms for \(N = 2, \dots, 60\), finds the lattice width by the same search over primitive directions the figure runs, and prints the bounds of the reformulated variable \(s\).
# Lenstra's strip K_N = {0 <= x <= N, 1 <= 10y - 5x <= 4}.
#
# The node counts of Proposition 1.5.16 for a branch and bound that
# branches on one integer form c.p (the search fig-lenstra runs), the
# lattice width by search, and the unimodular reformulation.
import math
def clip(poly, a, b, r):
"""Keep the part of a convex polygon with a x + b y <= r."""
out = []
for i in range(len(poly)):
P, Q = poly[i], poly[(i + 1) % len(poly)]
sP, sQ = a*P[0] + b*P[1] - r, a*Q[0] + b*Q[1] - r
if sP <= 1e-12:
out.append(P)
if (sP < -1e-12 and sQ > 1e-12) or (sQ < -1e-12 and sP > 1e-12):
t = sP / (sP - sQ)
out.append((P[0] + t*(Q[0] - P[0]), P[1] + t*(Q[1] - P[1])))
return out
def strip(N):
"""The vertices of K_N."""
return [(0, 0.1), (N, N/2 + 0.1), (N, N/2 + 0.4), (0, 0.4)]
def rng(poly, c):
"""The range (min, max) of the form c.p on a polygon."""
return (min(c[0]*p[0] + c[1]*p[1] for p in poly),
max(c[0]*p[0] + c[1]*p[1] for p in poly))
def holds_integer(lo, hi):
return math.ceil(lo - 1e-9) <= math.floor(hi + 1e-9)
def search(N, c):
"""Depth-first search along the form c.p; returns the node count.
A node is [k_lo, k_hi], with k_lo <= c.p <= k_hi.
"""
nodes, stack = 0, [(-math.inf, math.inf)]
while stack:
k_lo, k_hi = stack.pop()
nodes += 1
P = strip(N)
if k_hi < math.inf:
P = clip(P, c[0], c[1], k_hi)
if k_lo > -math.inf and P:
P = clip(P, -c[0], -c[1], -k_lo)
if len(P) < 1:
continue # relaxation empty
if (not holds_integer(*rng(P, (1, 0)))
or not holds_integer(*rng(P, (0, 1)))):
continue # closed by rounding x or y
lo, hi = rng(P, c)
if holds_integer(lo, hi):
# one child per lattice line c.p = k meeting the strip
ks = range(math.floor(hi + 1e-9), math.ceil(lo - 1e-9) - 1, -1)
stack.extend((k, k) for k in ks)
else:
# a single disjunction, both sides empty
stack.extend([(math.ceil(hi), math.inf),
(-math.inf, math.floor(lo))])
return nodes
forms = {"x": (1, 0), "y": (0, 1), "x + y": (1, 1), "2y - x": (-1, 2)}
print("N " + "".join(f"{name:>9}" for name in forms))
for N in (5, 10, 20, 40):
print(f"{N:<6}" + "".join(f"{search(N, c):>9}" for c in forms.values()))
closed = all(search(N, (1, 0)) == N + 2 and search(N, (-1, 2)) == 3
for N in range(2, 61))
print("closed forms N + 2 along x and 3 along 2y - x "
f"hold for N = 2..60: {closed}")
for N in (10, 40):
# the primitive directions (a, b) with |a|, |b| <= 6, one of each
# pair +-(a, b), by width of K_N along them, then by |a| + b
dirs = sorted({(a, b) for a in range(-6, 7) for b in range(0, 7)
if (a, b) != (0, 0) and math.gcd(abs(a), b) == 1
and (b > 0 or a > 0)},
key=lambda c: (round(rng(strip(N), c)[1]
- rng(strip(N), c)[0], 9),
abs(c[0]) + c[1]))
w = lambda c: rng(strip(N), c)[1] - rng(strip(N), c)[0]
print(f"N = {N}: lattice width {w(dirs[0]):.1f} along {dirs[0]},")
print(f" then {w(dirs[1]):.1f} along {dirs[1]}"
f" and {w(dirs[2]):.1f} along {dirs[2]}")
lo, hi = rng(strip(10), (-1, 2))
print("unimodular coordinates s = 2y - x, t = y: "
f"s in [{lo:.1f}, {hi:.1f}],")
print(f" rounded to [{math.ceil(lo)}, {math.floor(hi)}]: "
"empty before any branching")
Output:
N x y x + y 2y - x
5 7 3 8 3
10 12 6 16 3
20 22 11 31 3
40 42 21 61 3
closed forms N + 2 along x and 3 along 2y - x hold for N = 2..60: True
N = 10: lattice width 0.6 along (-1, 2),
then 5.3 along (0, 1) and 5.3 along (-1, 1)
N = 40: lattice width 0.6 along (-1, 2),
then 20.3 along (0, 1) and 20.3 along (-1, 1)
unimodular coordinates s = 2y - x, t = y: s in [0.2, 0.8],
rounded to [1, 0]: empty before any branching
(The cost per node, and three qualifications) Each node costs at most two polygon clippings and at most three range computations. The nodes of a search along a coordinate form are independent of each other once the root has been processed, so the whole tree could be bounded in one batch. The point of the example is that along the right form there is no tree to batch. Three remarks qualify the example. First, the strip is not close to the planar flatness constant. Its width \(0.6\) is below the slab bound \(1\) of Proposition 1.5.9, which is the right yardstick for a piece of a slab. The body that attains \(\mathrm{Flat}(2) = 1 + 2/\sqrt 3\) is Hurkens' triangle, which is wide in every direction, has a lattice point on each edge, and gives Lenstra's step two children rather than none. Second, the counts of Proposition 1.5.16 hold for \(N \ge 2\), and the figure's slider starts at \(5\). Third, branching on a general integer form \(\pi^\top x \le \pi_0\) or \(\pi^\top x \ge \pi_0 + 1\) is Lenstra's step with the width guarantee removed, and it is not a cure for exponential trees in varying dimension. Selecting the best such disjunction is itself NP-hard, and there are explicit families on which every branch-and-bound tree is exponential even when it may branch on arbitrary integer disjunctions.A. Mahajan and T. Ralphs, "On the complexity of selecting disjunctions in integer programming", SIAM Journal on Optimization 20 (2010); S. S. Dey, Y. Dubey and M. Molinaro, "Lower bounds on the size of general branch-and-bound trees", Mathematical Programming 198 (2023). Computational studies of the rule: J. H. Owen and S. Mehrotra, "Experimental results on using general disjunctions in branch-and-bound for general-integer linear programs", Computational Optimization and Applications 20 (2001); M. Karamanov and G. Cornuéjols, "Branching on general disjunctions", Mathematical Programming 128 (2011). What the theorem removes is the dependence on the data for a fixed number of integer variables, and that is all it removes.
A fixed number of integer variables, any number of continuous ones
Nothing in Theorem 1.5.14 used linearity. The argument needs a convex body, for the flatness, a rounding, the fact that a hyperplane section of a convex set is convex, for the recursion, and a bound on the region to be searched. All of this survives when the constraints are convex polynomials, and a little more. A function is quasiconvex if every sublevel set \(\{x : f(x) \le \alpha\}\) is convex. Every convex function is quasiconvex, and the converse fails, as the remark after Proposition 1.2.5 noted with \(\sqrt{|x|}\).
Theorem 1.5.17 (convex integer minimization in fixed dimension; Khachiyan and Porkolab 2000; Heinz 2005; Hildebrand and Köppe 2013). Let \(f, g_1, \dots, g_m\) be quasiconvex polynomials with integer coefficients, of degree at most \(d \ge 2\) and coefficient bit size at most \(\ell\), in \(n\) variables. There is an algorithm that computes a minimizer of \(\min\{ f(x) : g_i(x) \le 0,\ x \in \mathbb{Z}^n \}\), or reports that none exists, in time \(m\, \ell^{O(1)} d^{O(n)} 2^{O(n^3)}\), and the minimizer has bit size \(\ell\, d^{O(n)}\). Hildebrand and Köppe improve the dependence on \(n\) to \(2^{O(n \log n)}\), and Khachiyan and Porkolab prove polynomiality in fixed dimension for the more general problem of minimizing a convex polynomial over the integer points of a convex semialgebraic set.L. Khachiyan and L. Porkolab, "Integer optimization on convex semialgebraic sets", Discrete & Computational Geometry 23 (2000); S. Heinz, "Complexity of integer quasiconvex polynomial optimization", Journal of Complexity 21 (2005), with the running time as stated in Hemmecke, Köppe, Lee and Weismantel (2010), Theorem 10; R. Hildebrand and M. Köppe, "A new Lenstra-type algorithm for quasiconvex polynomial integer minimization with complexity \(2^{O(n \log n)}\)", Discrete Optimization 10 (2013). The oracle version of Lenstra's algorithm for convex bodies is in M. Grötschel, L. Lovász and A. Schrijver, Geometric Algorithms and Combinatorial Optimization (Springer, 1988).
(Convexity plus fixed dimension) Together with Theorem 1.5.11 this says that convexity plus fixed dimension is tractable, which is the integer analogue of the continuous fact that convex programming is polynomial. What changes between the three papers is where the rounding and the size bounds come from: real algebraic geometry for Khachiyan and Porkolab, direct estimates for Heinz, optimal ellipsoid roundings for Hildebrand and Köppe. The theorem fixes the total dimension. For the mixed-integer convex problems of this series, with a few integer variables and many continuous ones, the parameter one would like to count is the number of integers. The theory for that count is partial, and Basu's survey is the place to start.A. Basu, "Complexity of optimizing over the integers", arXiv 2110.06172 (2021). The partial results: T. Oertel, C. Wagner and R. Weismantel, "Integer convex minimization by mixed integer linear optimization", Operations Research Letters 42 (2014), minimize a convex function over integer points through a mixed-integer linear oracle; M. Baes, T. Oertel, C. Wagner and R. Weismantel, "Mirror-descent methods in mixed-integer convex optimization", in Facets of Combinatorial Optimization (Springer, 2013), adapt mirror descent to the mixed-integer setting; M. Baes, T. Oertel and R. Weismantel, "Duality for mixed-integer convex minimization", Mathematical Programming 158 (2016), develop a duality theory; A. Basu, H. Jiang, P. Kerger and M. Molinaro, "Information complexity of mixed-integer convex optimization", arXiv 2308.11153 (2023), bound the number of oracle calls any algorithm needs. For the linear case Lenstra's theorem makes the count exact: a mixed-integer linear program with a fixed number of integer variables is polynomial whatever the number of continuous ones.
(Two cautions) Two cautions belong here, because the result is often read as advice. The constants of these theorems are exponential in \(n\), \(2^{O(n \log n)}\) at best for the convex case, and no production solver implements a Lenstra-type recursion. The per-node work of a rounding, a reduction and several LPs buys a guarantee on the fan-out that variable branching with strong bounding does not need on typical instances. What solvers take from the theory is the thin-direction idea in the two forms just described, general-disjunction branching and lattice reformulation as preprocessing. And the theorems are about convex sets. The next part shows that none of this survives without convexity.
Nonconvexity is not rescued by fixed dimension
A fully polynomial-time approximation scheme is an algorithm that, for every \(\varepsilon > 0\), returns a solution within a factor \(1 + \varepsilon\) of the optimum in time polynomial in the input size and in \(1/\varepsilon\).
Theorem 1.5.18 (polynomial objectives; De Loera, Hemmecke, Köppe and Weismantel 2006, 2008; Hemmecke, Köppe, Lee and Weismantel 2010). (i) Minimizing a polynomial of degree \(4\) over the integer points of a convex polygon, that is, in dimension \(2\), is NP-hard. (ii) Continuous polynomial optimization over polytopes in varying dimension is NP-hard, and admits no fully polynomial-time approximation scheme unless P = NP. (iii) For fixed dimension there is a fully polynomial-time approximation scheme for maximizing a polynomial that is nonnegative over the integer points of a polytope. The same holds for the mixed-integer points of a polytope when the number of continuous variables is also fixed.J. A. De Loera, R. Hemmecke, M. Köppe and R. Weismantel, "Integer polynomial optimization in fixed dimension", Mathematics of Operations Research 31 (2006); J. A. De Loera, R. Hemmecke, M. Köppe and R. Weismantel, "FPTAS for optimizing polynomials over the mixed-integer points of polytopes in fixed dimension", Mathematical Programming 115 (2008); statements (i) and (ii) in the form given are Theorems 2 and 1 of Hemmecke, Köppe, Lee and Weismantel (2010). For quadratics: R. Hildebrand, R. Weismantel and K. Zemmer, "An FPTAS for minimizing indefinite quadratic forms over integers in polyhedra", SODA 2016, and, most recently, A. Del Pia, "Rational Jacobi rotations and the complexity of approximating mixed integer quadratic programming", arXiv 2607.29386 (2026).
(Where flatness breaks without convexity) Part (i) is the sharpest possible contrast with Theorem 1.5.17: two integer variables, a polygon and degree four already encode an NP-complete problem, the reduction being from a quadratic congruence problem. Part (ii) is max-cut, the problem of splitting a graph's vertices into two sets so as to maximize the number of edges between them, written as \(\min x^\top Q x\) over \([-1, 1]^n\). Behind the second clause is Håstad's inapproximability theorem, stated in Theorem 1.5.19 below. Part (iii) says that nonnegativity plus approximation recovers something in fixed dimension, through the evaluation of power sums over the lattice points of the polytope rather than through flatness. Where the flatness argument breaks is easy to see. The flatness theorem is a statement about convex sets, and a lattice-point-free set that is not convex need not be thin in any direction. The set \(\{ x \in [0, N]^2 : \|x - m\|_\infty \ge \tfrac14 \text{ for all } m \in \mathbb{Z}^2 \}\), a square with a small hole punched at every lattice point, contains no lattice point and has width about \(N\) in every direction. The feasible set of a single quadratic equation \(w = y_1 y_2\), which is the gadget behind Theorem 1.5.4, is of this kind, and so is the sublevel set of the quartic in (i). The continuous side behaves differently: a nonconvex polynomial program in fixed dimension is polynomial, because Renegar's \((sd)^{O(n)}\) is a polynomial for fixed \(n\), so the hardness in (i) is a joint effect of the integers and the curve, neither alone. Quadratic objectives sit on the boundary. Mixed-integer quadratic programming is in NP (Theorem 1.5.7), the two-variable integer case is polynomial, and polynomiality of integer quadratic programming in fixed dimension is open.
Approximation does not rescue it either, and convexity is a declaration
A reader who has followed the nonconvex case this far might hope that if exact optimization is out of reach, a relaxation with a guaranteed ratio is not. For quadratic programs it is.
Theorem 1.5.19 (hardness of approximation; Bellare and Rogaway 1995; Håstad 2001). Maximizing a polynomial over a polytope admits no polynomial-time approximation algorithm with a good guarantee unless P = NP, where "good" means a relative error measured against the range of the objective over the polytope. Quadratic programming admits none unless NP is contained in quasi-polynomial time, that is, time \(2^{\mathrm{polylog}}\) in the input size. For max-cut, approximating the optimum within a factor \(16/17 + \varepsilon\) is NP-hard for every \(\varepsilon > 0\).M. Bellare and P. Rogaway, "The complexity of approximating a nonlinear program", Mathematical Programming 69 (1995), whose paper makes "good" precise; J. Håstad, "Some optimal inapproximability results", Journal of the ACM 48 (2001). The positive side for max-cut, the semidefinite relaxation with a \(0.878\) guarantee of Goemans and Williamson, is in Section 4.7.
(No relaxation is uniformly tight) The consequence for the design of relaxations, which is the subject of Sections 2 and 4, is that no relaxation of nonconvex quadratic programs can be uniformly tight. A convex relaxation whose value was always within a fixed factor of the optimum would estimate the optimal value within that factor, which the promise-problem versions of these theorems (gap problems in the complexity sense) forbid. Every relaxation in this series is therefore tight on some instances and loose on others, and the search of Section 3 exists to pay back the difference on the instances where it is loose.
(Convex is a declaration, not a discovery) A last fact closes the circle with Section 1.2. The watershed between convex and nonconvex problems is real, but a solver cannot locate it by inspection. Burer's theorem concerns nonconvex quadratic programs with linear equality constraints, nonnegative variables and binary variables, under his key assumption that the linear constraints bound the binary variables between \(0\) and \(1\). It says that every such program is equal in value to a linear program over the completely positive cone, the matrices that are sums of \(x x^\top\) with \(x \ge 0\). The nonconvexity has been moved entirely into a convex cone. Membership in the completely positive cone is NP-hard, and membership in its dual, the copositive cone of Theorem 1.3.9, is co-NP-complete. The reformulation is therefore exact and intractable, and "convex" in the sense of "written as a convex program" is not the same as "easy".S. Burer, "On the copositive representation of binary and continuous nonconvex quadratic programs", Mathematical Programming 120 (2009); P. J. C. Dickinson and L. Gijben, "On the computational complexity of membership problems for the completely positive cone and its dual", Computational Optimization and Applications 57 (2014). The tractable inner approximations of the cone are the subject of Section 4.7. In the other direction, Section 1.2 recorded that deciding the convexity of a quartic is NP-hard (Ahmadi, Olshevsky, Parrilo and Tsitsiklis 2013, stated in full as Theorem 2.5.15), while for quadratics it is a positive semidefiniteness test. A "convex MINLP" is therefore a syntactic class: a problem whose convexity is established by construction. The construction is a positive semidefiniteness test for quadratics, the composition rules of disciplined convex programming, a grammar of operations that preserve convexity, or the tree walks over the expression graph of Section 2.5. A nonconvex solver does not discover hidden convexity in general, and the convex methods of Section 3.4 are applied to problems that declare it.
The table
The four-letter table of Section 1.1 sorted problems by what the functions and variables may be. The following table sorts the same classes by what is known, with the theorems of this subsection as its entries. "Bounded" means finite bounds on every variable, and the data are rational throughout.
| class | decidable? | in NP? | NP-hard? | fixed dimension | exact rational optimum? |
|---|---|---|---|---|---|
| LP | yes | yes: in P (Khachiyan 1980) | no | in P | yes: a vertex (Cramer's rule) |
| MILP, bounded or not | yes | yes (Papadimitriou 1981): NP-complete | yes (Karp 1972) | polynomial (Lenstra 1983) | yes: a vertex of \(\operatorname{conv}(\mathcal F)\), the integer hull of 2.1 (Meyer 1974) |
| convex NLP, polynomial data | yes (Tarski; Renegar 1992) | open in general; convex QP is in P | no: an \(\varepsilon\)-optimum in polynomial time | polynomial to accuracy \(\varepsilon\) in any dimension (ellipsoid) | no: \(\min -x\) s.t. \(x^2 \le 2\) has optimum \(-\sqrt 2\) |
| nonconvex NLP, polynomial data | yes (Tarski; Renegar 1992) | feasibility is \(\exists\mathbb{R}\)-complete; NP membership open | yes, already with one negative eigenvalue (Pardalos–Vavasis 1991) | polynomial for fixed \(n\) (Renegar: \((sd)^{O(n)}\)) | no |
| MIQP: nonconvex quadratic objective, linear constraints | yes | yes (Del Pia, Dey and Molinaro 2017) | yes | open; polynomial for two integer variables | yes (bounded): a KKT point of a rational linear system |
| convex MINLP, bounded | yes (enumerate the integers; Renegar) | open | yes (contains MILP) | polynomial (Khachiyan–Porkolab 2000; Heinz 2005) | no |
| nonconvex MINLP, bounded, polynomial data | yes (enumerate the integers; Renegar) | open | yes | NP-hard already for degree 4 in the plane (De Loera et al.) | no |
| MINLP, unbounded integers, quadratic constraints | NO (Jeroslow 1973) | – | – | undecidable with 10 integer variables (Jones 1982) | – |
| NLP or MINLP with exp, log, sin, unbounded variables | not known (Richardson 1968) | – | – | – | – |
(Three readings of the table) Three readings of the table matter for the rest of the series. Down the decidability column, the only thing that destroys decidability is unbounded integers together with a product of variables. Once every variable is bounded, every polynomial problem on the list is decidable, so the bounds the solvers insist on are what keep them in the decidable rows. Down the fixed-dimension column, fixed dimension helps exactly where there is convexity, and the plane with a quartic is already as hard as anything. Down the last column, an exact rational answer exists only where the constraints are linear, in the LP, MILP and MIQP rows. For a quadratic objective over a polyhedron the minimum is attained at a KKT point. The KKT system of its active set is a rational linear system, and the objective is constant on that system's solution set, so a rational optimal point exists. With bounded integer variables the same holds for MIQP. Once a constraint is nonlinear the optimum is algebraic in general, which is why every nonlinear solver works to a tolerance. The exact-arithmetic solvers of Section 5.6 exist today only for MILP. For MIQP that is a gap in the tooling rather than an impossibility. The entries of the table not attributed above are the polynomiality of linear programming and of convex quadratic programming, and the ellipsoid method behind the convex row.L. G. Khachiyan, "Polynomial algorithms in linear programming", USSR Computational Mathematics and Mathematical Physics 20 (1980); N. Karmarkar, "A new polynomial-time algorithm for linear programming", Combinatorica 4 (1984); M. K. Kozlov, S. P. Tarasov and L. G. Khachiyan, "The polynomial solvability of convex quadratic programming", USSR Computational Mathematics and Mathematical Physics 20 (1980); M. Grötschel, L. Lovász and A. Schrijver, "The ellipsoid method and its consequences in combinatorial optimization", Combinatorica 1 (1981); M. Frank and P. Wolfe, "An algorithm for quadratic programming", Naval Research Logistics Quarterly 3 (1956), for the attainment of the minimum of a quadratic bounded below on a polyhedron.
What a solver certifies
That last point is the one that makes solvers terminate. The optimum of \(\min -x\) subject to \(x^2 \le 2\) is \(-\sqrt 2\), a convex problem with rational data whose answer no finite decimal writes down, so a solver that promised the exact optimum could not stop. By contrast \(\min x\) subject to \(3x \ge 2\) has the rational answer \(\tfrac23\). What a global solver returns is the following.
Definition 1.5.20 (what a global solver certifies). Given tolerances \(\varepsilon_{\mathrm{feas}}, \varepsilon_{\mathrm{int}}, \varepsilon_a, \varepsilon_r > 0\), a global solver returns a point \(\bar x\) and a number \(\underline z\) such that (i) \(g_i(\bar x) \le \varepsilon_{\mathrm{feas}}\) for every \(i\); (ii) \(|\bar x_j - \operatorname{round}(\bar x_j)| \le \varepsilon_{\mathrm{int}}\) for every integer coordinate \(j\); (iii) \(\underline z \le z^\star\) is a valid lower bound; and (iv) \(f(\bar x) - \underline z \le \max(\varepsilon_a,\ \varepsilon_r |f(\bar x)|)\). The run is a certificate that a nearly feasible point is \(\varepsilon\)-globally optimal.
On R3 of Section 1.6, for instance, the run of Section 3.5 ends with \(\bar x = (0.5, 1)\), which satisfies both constraints exactly and has an integral \(y\), so (i) and (ii) hold with zero slack, and with a lower bound \(\underline z\) within the run's gap tolerance of \(f(\bar x) = 2.0\), so (iii) and (iv) hold. On most instances the point returned is only nearly feasible, its integer coordinates are only nearly integral, and the gap stops at the tolerance rather than at zero, which is why all four tolerances appear in the definition.
(Why the tolerances are part of the problem) Condition (iii) is the one every bound computed in a tree has to respect, including in floating-point arithmetic, and Sections 7.3 and 8 return to it. The existence of such a certificate for every bounded MINLP with continuous data is the convergence theorem for spatial branch and bound of Section 3.5. With an exhaustive subdivision of the box and consistent bounds, the gap between the incumbent and the lowest open bound tends to zero, and it falls below any positive tolerance after finitely many nodes.R. Horst and H. Tuy, Global Optimization: Deterministic Approaches, 3rd ed. (Springer, 1996), Chapter IV, "Branch and Bound", pp. 115–178. Section 2.1 tabulates the gap conventions and the default gap tolerances, and Section 8 returns to all three tolerances. None of the exact decision problems of this subsection is ever posed by a solver, and the tolerances are not an engineering compromise laid over an exact method. They are part of the problem statement.
Where this is used
Every global solver in Section 5 requires finite bounds on the variables that appear in nonconvex terms, and infers them when the model does not provide them. It uses them for far more than termination, because every envelope of Section 2.4 and every propagation of Section 2.6 is built on a box. That is Proposition 1.5.6 and Theorem 1.5.7 in production. No solver implements Lenstra's recursion, but the two ideas it leaves behind are in use: lattice reformulation for equality-constrained integer problems, and general-disjunction branching as an option. The three kinds of tolerance of Definition 1.5.20, feasibility, integrality and gap (absolute or relative), are the ones every solver's log reports, under names Section 8 tabulates.
What parallelizes
For parallel and GPU work the results of this subsection draw one line. On the hard families of Section 3.1, the work of a branch-and-bound run is exponential in the number of integer variables whatever the order in which nodes are processed. A machine that processes nodes a thousand times faster therefore leaves the exponent unchanged: since \(2^{10} = 1{,}024\), it moves the largest instance that fits a fixed budget by about ten binary variables. What changes the base of the exponential is the relaxation, the propagation, the cuts and the disjunctions branched on, which is why the rest of this series is about those. Lenstra's theorem separates the two sources of difficulty cleanly. The data enter only through polynomial-time work at each node, which is where parallel hardware applies, while the dimension enters through the fan-out of the tree, which it does not touch. The programme of Section 7 is to make the per-node work, the bound, cheap and massively parallel, in the knowledge that the number of nodes is set by the mathematics of the previous pages and not by the hardware.
The running examples
Six examples recur through the series, and the preceding subsections have already pointed to several of them. Each is small enough to compute by hand or with a few lines of code, and each is chosen to isolate one phenomenon. This subsection gives their complete data, the figures in which they appear, and the headline numbers the later sections quote, so that any number in the text can be traced back to one of these six definitions. It closes this section because Sections 2 to 4 quote these numbers without re-deriving them. Two programs at the end reproduce the central numbers of the six examples. The remaining numbers are quoted from the live figures and from the computations named in the sidenotes.
R1. A two-variable integer program over a polygon. The drawn example of the plane figures is
\[\max_{x,\,y \in \mathbb{Z}}\ c_1 x + c_2 y \quad\text{subject to}\quad 2x + 5y \le 24.5,\ \ 5x + 2y \le 30.5,\ \ -3x + 4y \le 11,\ \ x - 2y \le 4.2,\ \ 0 \le x \le 7,\ \ 0 \le y \le 6,\]with the objective direction \(c = (\cos\theta, \sin\theta)\) set by an angle the figures let the reader turn. It maximizes, and the figures say so. In the notation of the displays it minimizes \(-c^\top (x, y)\). The polygon \(P\) cut out by the four rows and the nonnegativity constraints has area \(18.245\) and contains \(22\) integer points, whose convex hull \(H\) has \(7\) facets. At the default angle \(\theta = 45°\) the LP relaxation over \(P\) reaches \(5.556\) at the fractional vertex \((4.929, 2.929)\). A best integer point is \((4, 3)\), tied with \((5, 2)\), with value \(4.950\), so the relaxation overstates the optimum by \(0.606\), which is \(10.9\%\) of the bound. A second polygon \(W\), each facet of \(P\) pushed outward until the next integer point would enter, describes the same 22 points with a weaker relaxation: \(5.627\) at \((4.979, 2.979)\) and a gap of \(12.0\%\), the numbers the relaxation figure's readout prints. The example is for everything that can be drawn in a plane. Section 2.1 uses it for the relaxation and the gap, and Sections 2.1 and 2.3 for the integer hull and the vocabulary of a solve. Section 2.6 tightens its box \([0, 7] \times [0, 6]\). Section 3.3 adds cutting planes at its root, and Section 3.2 runs the feasibility pump and the neighbourhood searches on it. Section 7.2 solves its LP relaxation three ways, by the simplex method, by a barrier method and by PDHG. Proposition 1.4.1 and the paragraph after it show that its LP relaxation is the convex-envelope relaxation of the two integrality constraints, written there as \(g(x) \le 0\) and \(g(y) \le 0\), so replacing each by its envelope leaves exactly \(P\).
R2. The bilinear example. The smallest problem with a product of two variables is
\[\max\ xy \quad\text{subject to}\quad 2x + y \le 1.2, \qquad (x, y) \in [0, 1]^2 .\]It maximizes. The constraint is active at the optimum, so the problem is \(\max\{x(1.2 - 2x) : 0 \le x \le 0.6\}\), whose derivative \(1.2 - 4x\) vanishes at \(x = 0.3\). The answer is \(0.18\) at \((0.3, 0.6)\), where the line is tangent to the level curve \(xy = 0.18\) and the gradient \((0.6, 0.3)\) of the product is parallel to the normal \((2, 1)\) of the line. The McCormick relaxation of Section 2.4, which replaces the product \(xy\) by a new variable held between four planes, answers \(0.400\) at \((0.4, 0.4)\), where the true product is only \(0.16\). Section 2.4 derives the four planes and this value.G. P. McCormick, "Computability of global solutions to factorable nonconvex programs: Part I—Convex underestimating problems", Mathematical Programming 10 (1976), for the relaxation; F. A. Al-Khayyal and J. E. Falk, "Jointly constrained biconvex programming", Mathematics of Operations Research 8 (1983), for the proof that it is the tightest possible on a rectangle. The example is for relaxation strength and for what closes a gap. Partitioning each axis into \(k\) pieces and relaxing each sub-box separately brings the relaxed maximum down to \(0.233\), \(0.193\), \(0.192\), \(0.187\) and \(0.185\) for \(k = 2, \dots, 6\), at the price of \(k^2\) boxes (Section 2.4). Spatial branch and bound on the same relaxation proves \(0.18\) to a tolerance of \(0.01\) in \(13\) nodes with midpoint splits, to \(0.001\) in \(21\) nodes and to \(0.0001\) in \(27\). A uniform grid needs \(k = 5\), that is \(25\) boxes, for the first tolerance (Section 3.5). Multiplying the constraints by one another and by the bounds, the reformulation-linearization technique of Section 4.7, lowers the bound at the root, the first relaxation a search solves before any branching, from \(0.400\) to \(0.300\). One further convexity fact about the square of \(x\), stated in Section 4.7, lowers it to \(0.180\), the exact value, with no branching at all.H. D. Sherali and A. Alameddine, "A new reformulation-linearization technique for bilinear programming problems", Journal of Global Optimization 2 (1992). The figure of Section 4.7 solves the third linear program with tangents of \(x^2\) at nine points and reports the vertex \((x, y, W, X, Y) = (0.25, 0.55, 0.18, 0.06, 0.30)\) of a flat optimal face on which \((0.3, 0.6, 0.18, 0.09, 0.36)\) also lies; the value is \(0.180\) either way.
| way | effort | bound on max xy |
|---|---|---|
| McCormick relaxation of the whole square (k = 1) | 1 box | 0.400 |
| uniform grid, k = 2 pieces per axis | 4 boxes | 0.233 |
| uniform grid, k = 3 | 9 boxes | 0.193 |
| uniform grid, k = 4 | 16 boxes | 0.192 |
| uniform grid, k = 5 | 25 boxes | 0.187, within 0.01 |
| uniform grid, k = 6 | 36 boxes | 0.185 |
| spatial branch and bound, midpoint splits | 13 nodes | within 0.01 |
| spatial branch and bound, midpoint splits | 21 nodes | within 0.001 |
| spatial branch and bound, midpoint splits | 27 nodes | within 0.0001 |
| RLT at the root | no branching | 0.300 |
| RLT and the convexity of x² at the root | no branching | 0.180 |
R3. MINLPLib st_e13. A two-variable library instance on which a convex-case method fails is
\[\min\ 2x + y \quad\text{subject to}\quad 1.25 - x^2 - y \le 0, \qquad x + y \le 1.6, \qquad 0 \le x \le 1.6, \qquad y \in \{0, 1\}.\]The nonlinear constraint reads \(y \ge 1.25 - x^2\). Its feasible side is the region on or above a downward parabola, so the constraint function is concave and the feasible set is two horizontal segments: \(y = 1\) with \(x \in [0.5, 0.6]\), and \(y = 0\) with \(x \in [\sqrt{1.25}, 1.6]\). On each segment the objective is smallest at the left end. The two candidates are therefore \((0.5, 1)\) with value \(2.0\), the optimum, and \((1.118, 0)\) with value \(2\sqrt{1.25} = \sqrt 5 = 2.236\). MINLPLib records the instance as type MBQCP, a mixed-binary quadratically constrained program, with two variables (one binary) and two constraints (one linear, one quadratic, curvature concave). It lists the primal bound \(2\) and the dual bound \(2\), the latter reported by ANTIGONE, BARON, Couenne, Gurobi, LINDO and SCIP, names the source "BARON book instance misc/e13", and dates the instance's addition to 1 September 2002.MINLPLib, instance page st_e13, minlplib.org/st_e13.html, read 5 October 2026. The library names the source as the problem directory misc/e13 of M. Tawarmalani and N. V. Sahinidis, Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming (Kluwer, 2002), and lists as a second reference G. R. Kocis and I. E. Grossmann, "Global optimization of nonconvex mixed-integer nonlinear programming (MINLP) problems in process synthesis", Industrial & Engineering Chemistry Research 27 (1988); whether the instance appears in the latter was not checked. The GAMS model on the page writes the quadratic row as \(-x^2 - y \le -1.25\). Three results from it recur. Outer approximation, the method of Section 3.4 for convex problems, started from \(y = 0\) finds the point \((1.118, 0)\) and linearizes the parabola there into the cut \(y \ge 2.5 - 2.236\,x\). The optimum violates that cut, because a tangent of a concave function lies above it (Proposition 1.2.4). The method then declares the subproblem with \(y = 1\) infeasible and stops at \(2.236\) with a gap of zero and a wrong answer.M. A. Duran and I. E. Grossmann, "An outer-approximation algorithm for a class of mixed-integer nonlinear programs", Mathematical Programming 36 (1986); R. Fletcher and S. Leyffer, "Solving mixed integer nonlinear programs by outer approximation", Mathematical Programming 66 (1994), whose finiteness theorem assumes convexity. At the optimum the cut reads \(2.236 \cdot 0.5 + 1 = 2.118 < 2.5\). Replacing \(-x^2\) on \([0, 1.6]\) by its chord gives the relaxed constraint \(y \ge 1.25 - 1.6x\) and the root bound \(1.3125\) at \((0.15625, 1)\). Fixing \(y = 1\) leaves that bound at \(1.3125\), so branching on the integer variable alone cannot close the gap. The integrality gap of an instance is what remains of the gap when only the integrality of \(y\) is dropped, and its nonconvexity gap is what remains when the curved constraint is relaxed as well. Here the first is zero and the whole gap is a nonconvexity gap. Spatial branch and bound with chord relaxations, one round of bound tightening, the inference of tighter variable bounds from the constraints on the current box (Section 2.6), before each spatial branch and a split of \(x\) at the incumbent's value proves \(2.0\) in \(5\) nodes, \(7\) relaxations and \(2\) tightening rounds. These are the counts the figure of Section 3.5 displays. A hand trace of the same search that skips the relaxations whose outcome the incumbent already decides counts five relaxations. Both are correct runs of the same algorithm under two counting conventions. The counts describe the figure's deliberately simple tightening rule. The figure's other two tightening settings give \(7\) nodes and \(7\) relaxations with no tightening and \(3\) nodes and \(9\) relaxations with tightening run to a fixed point (Section 3.5). A solver propagates the quadratic constraint itself: once \(y\) is fixed it reads \(x \ge \sqrt{1.25 - y}\) off the constraint directly, as Couenne, SCIP and BARON do as a matter of course. That step is not one of the figure's settings. A variant of the figure's algorithm with it, run for this series by the script of Section 3.5, closes each node with \(y\) fixed by its own relaxation and proves the optimum in \(3\) nodes and \(3\) relaxations.Feasibility-based tightening through the nonlinear constraints is part of each solver's node loop: H. S. Ryoo and N. V. Sahinidis, "A branch-and-reduce approach to global optimization", Journal of Global Optimization 8 (1996), for BARON's lineage; Belotti, Lee, Liberti, Margot and Wächter (2009), cited in Section 1.3, for Couenne; K. Bestuzheva, A. Chmiela, B. Müller, F. Serrano, S. Vigerske and F. Wegscheider, "Global optimization of mixed-integer nonlinear programs with SCIP 8", Journal of Global Optimization 91 (2025), for SCIP. Section 2.6 gives the algorithm.
R4. Haverly's pooling problem. Three inputs with different sulphur contents are blended for two products with sulphur limits. Inputs A and B enter a pool, input C bypasses it, and the pool sends its mixture to both products:
| input | cost | sulphur % | path | product | price | sulphur limit % | demand (at most) |
|---|---|---|---|---|---|---|---|
| A | 6 | 3 | into the pool | X | 9 | 2.5 | 100 |
| B | 16 | 1 | into the pool | Y | 15 | 1.5 | 200 |
| C | 10 | 2 | direct to X, Y |
With flows \(f_A, f_B\) into the pool, \(y_X, y_Y\) from the pool to the products, \(z_X, z_Y\) from C to the products, and the pool's sulphur content \(q \in [1, 3]\), the problem is
\[\begin{aligned} \max\;& 9\,(y_X + z_X) + 15\,(y_Y + z_Y) - 6 f_A - 16 f_B - 10\,(z_X + z_Y) \\ \text{s.t.}\;& f_A + f_B = y_X + y_Y, \qquad 3 f_A + f_B = q\,y_X + q\,y_Y, \\ & q\,y_X + 2 z_X \le 2.5\,(y_X + z_X), \qquad q\,y_Y + 2 z_Y \le 1.5\,(y_Y + z_Y), \\ & y_X + z_X \le 100, \qquad y_Y + z_Y \le 200, \qquad \text{all flows} \ge 0 . \end{aligned}\]Haverly's pooling network, instance 1
f_A y_X
A 6, 3 ---------+ +-------------------> X 9, 2.5, 100
| | ^
v | | z_X
+-----------------------+ |
| pool: sulphur q, | C 10, 2
| q in [1, 3] | |
+-----------------------+ | z_Y
^ | v
B 16, 1 ---------+ +-------------------> Y 15, 1.5, 200
f_B y_Y
inputs A, B, C: cost and sulphur %; products X, Y: price,
sulphur limit % and demand (at most)
at the pool f_A + f_B = y_X + y_Y
3 f_A + f_B = q y_X + q y_Y the two bilinear terms
It maximizes. The two bilinear terms \(q\,y_X\) and \(q\,y_Y\) are its only nonconvexity, and the pool quality \(q\) is a continuous variable, so there is nothing to branch on but a continuous variable.C. A. Haverly, "Studies of the behavior of recursion for the pooling problem", ACM SIGMAP Bulletin 25 (1978). The instance is MINLPLib's haverly, with optimum \(-400\) in minimization form. The instance-2 and instance-3 modifications are the ones used in the literature and reproduce the known optima \(600\) and \(750\); they were not checked against Haverly's original tables. The formulation above is the p-formulation; the q- and pq-formulations in proportions are Section 3.5's. For a fixed \(q\) the problem is linear in the flows, and sweeping \(q\) draws a landscape with two hills. The profit is \(400\) at \(q = 1\), where the pool takes only B and feeds Y, and C also feeds Y. It falls to \(300\) at \(q = 1.5\) and drops to \(0\) just above \(1.5\), where product Y can no longer be blended to its sulphur limit. It stays at \(0\) up to \(q = 2.4\), where the pool's material costs exactly what product X pays for it, then rises to \(50\) at \(q = 2.5\) and to \(100\) at \(q = 3\), where the pool takes only A and feeds X. The curve is a maximum over linear programs whose feasible sets change with \(q\), and the drop at \(q = 1.5\) is a genuine discontinuity. The global optimum is \(400\). The other hill, \(100\), is a local optimum, and Haverly's original method, which fixes \(q\), solves the LP and re-estimates \(q\) from the flows, stays there if it starts there. The McCormick relaxation of the two products claims \(500\), at \(q = 2\). Section 2.4 reads off its relaxed point why the bound is loose, and Section 3.5 draws it. The pq-formulation, which adds the products of the pool's proportion constraint with the outflows, gives the same bound \(500\) on this instance, although its relaxation is at least as tight as the p-formulation's in general. Splitting \(q\) once at \(2\%\) and relaxing each half gives \(400\) on \([1, 2]\) and \(100\) on \([2, 3]\), so one spatial branch on the pool quality proves the optimum (Sections 2.4 and 3.5). The order-2 moment relaxation of Section 4.7 closes the gap at the root instead, at \(400.0003\), with a \(36 \times 36\) matrix variable.The pq-formulation and the proof that its relaxation is at least as tight as the p-formulation's are in Tawarmalani and Sahinidis (2002), cited above, chapter 9; the moment relaxation is J. B. Lasserre, "Global optimization with polynomials and the problem of moments", SIAM Journal on Optimization 11 (2001). The bounds \(500\), \(1{,}000\) and \(800\) for the three instances were recomputed for this series with HiGHS 1.15.1 and again with a numpy-only simplex; the value \(400.0003\) was computed with cvxpy 1.9.3 and the conic interior-point solver Clarabel 0.11.1, with the variables scaled to \([0, 1]\). Section 3.5's figure computes the first bound and the split bounds live. Instances 2 and 3 move the optimum to the other hill and to the interior. Instance 2 has optimum \(600\) at \(q = 3\), a local optimum of \(400\) at \(q = 1\) and a root bound of \(1{,}000\). Instance 3 has optimum \(750\) at \(q = 1.5\), a local optimum of \(125\) at \(q = 2.5\) and a root bound of \(800\). The example is for the bilinear nonconvexity in its original setting. It shows local optima that a local method cannot see past, and spatial branching on a continuous variable (Section 3.5).
haverly; the instance-2 and instance-3 modifications were not checked against Haverly's original tables, and the bounds were recomputed with HiGHS 1.15.1 and a numpy-only simplex, as the sidenote says.R5. Three assets, at most two held. The quadratic running example is a mean–variance portfolio with a limit on the number of positions. Three assets have expected returns \(\mu = (10\%, 7\%, 4\%)\), volatilities \(\sigma = (20\%, 15\%, 10\%)\) and a common correlation \(\rho = 0.3\) between every pair, so the covariance is the one-factor matrix \(\Sigma = \rho\, \sigma\sigma^\top + (1 - \rho)\operatorname{diag}(\sigma^2)\). The weights \(w\) are long only and sum to one, the return must reach a target \(R\), and at most two assets may be held:
\[\min_{w,\,z}\ w^\top \Sigma w \quad\text{subject to}\quad e^\top w = 1, \qquad \mu^\top w \ge R, \qquad 0 \le w_i \le z_i, \qquad e^\top z \le 2, \qquad z \in \{0, 1\}^3 .\]The continuous relaxation of the indicator constraints, \(z \in [0, 1]^3\), allows \(z = w\) for every \(w\) on the simplex, since then \(0 \le w_i \le z_i\) holds and \(e^\top z = e^\top w = 1 \le 2\). The relaxation therefore forgets the limit entirely and returns the unconstrained minimum-variance portfolio. At the default target \(R = 7\%\) that portfolio has volatility \(11.12\%\). The best portfolio of two assets has \(12.45\%\) and holds assets 1 and 3 in equal halves. The perspective relaxation of Section 4.3, which treats a diagonal part of \(\Sigma\) exactly, gives \(11.71\%\) and closes \(43\%\) of the gap measured on the variances.The cardinality-constrained mean–variance problem is D. Bienstock, "Computational study of a family of mixed-integer quadratic programming problems", Mathematical Programming 74 (1996); the mean–variance problem itself is H. Markowitz, "Portfolio selection", The Journal of Finance 7 (1952). The figure of Section 4.4 extracts the diagonal \(d\,I\) with \(d = 0.99\,\lambda_{\min}(\Sigma)\); in variances the three numbers are \(0.012375\), \(0.015500\) and \(0.0137143\), and \((0.0137143 - 0.012375)/(0.015500 - 0.012375) = 0.429\). The example is for the indicator structure that pervades portfolio and tax problems. It shows the difference between a big-M formulation and a perspective formulation of the same logical constraint (Sections 4.2 to 4.4), and the way a factor model hands the perspective the diagonal it needs (Sections 4.4 and 4.5). It is also the smallest instance of the problem that Section 4.9 scales to thousands of securities.
| at R = 7% | volatility | variance |
|---|---|---|
| continuous relaxation: the limit forgotten | 11.12% | 0.012375 |
| perspective relaxation | 11.71% | 0.013714 |
| optimum with at most two assets: assets 1 and 3 in equal halves | 12.45% | 0.015500 |
R6. Eight tax lots. An account holds eight lots of one stock, bought on different dates at different prices, and must raise \(C = \$62{,}000\) by selling whole lots at today's price of \(\$100\) a share. Selling lot \(\ell\) with \(n_\ell\) shares and basis \(\beta_\ell\) realizes a gain of \(n_\ell\,(100 - \beta_\ell)\). The gain is taxed at \(40.8\%\) if the lot has been held for a year or less and at \(23.8\%\) if longer, the top federal rates plus the \(3.8\%\) net investment income tax. A loss carries a negative tax, which is a benefit.
| lot | shares | basis | days held | term | gain | tax |
|---|---|---|---|---|---|---|
| 1 | 100 | $62 | 900 | LT | 3,800 | 904.40 |
| 2 | 150 | $118 | 400 | LT | −2,700 | −642.60 |
| 3 | 80 | $95 | 340 | ST | 400 | 163.20 |
| 4 | 120 | $131 | 45 | ST | −3,720 | −1,517.76 |
| 5 | 60 | $88 | 500 | LT | 720 | 171.36 |
| 6 | 200 | $104 | 200 | ST | −800 | −326.40 |
| 7 | 90 | $71 | 1,200 | LT | 2,610 | 621.18 |
| 8 | 50 | $99 | 30 | ST | 50 | 20.40 |
With \(t_\ell\) the tax of lot \(\ell\) and \(p_\ell = 100\, n_\ell\) its proceeds, the problem is a knapsack problem: choose items of given weights to meet a capacity at least cost. Here the items are the lots, the weights are their proceeds, the capacity is the cash and the cost is the tax,
\[\min_{x \in \{0, 1\}^8}\ \sum_{\ell = 1}^{8} t_\ell\, x_\ell \quad\text{subject to}\quad \sum_{\ell = 1}^{8} p_\ell\, x_\ell \;\ge\; C .\]Of the \(256\) subsets of lots, \(34\) raise enough cash. The best of them sells lots 2, 3, 4, 5, 6 and 8, raises \(\$66{,}000\) and pays \(-\$2{,}131.80\) in tax, a net loss harvested. The three loss lots raise \(\$47{,}000\) and are worth \(\$2{,}486.76\). The cheapest gains per dollar of proceeds are lots 8, 3 and 5 in that order, and taking all three raises \(\$66{,}000\) against the \(\$62{,}000\) required. The two rules of thumb, selling the lots with the least tax per dollar of proceeds first and selling the highest-basis lots first, happen to choose the same six lots here. The LP relaxation replaces each \(x_\ell \in \{0, 1\}\) by its convex hull \(0 \le x_\ell \le 1\). Equivalently, by Proposition 1.4.1, it replaces the constraint \(x_\ell - x_\ell^2 \le 0\) by its convex envelope on \([0, 1]\), the zero function. It sells the loss lots whole, then lots 8 and 3, then a third of lot 5, raises exactly \(\$62{,}000\) and promises \(-\$2{,}246.04\). The \(\$114.24\) between the two is the integrality of lots, the price of selling all of lot 5 rather than a third of it, and it is the gap a solver would have to close. Move the calendar forward thirty days and lot 3 crosses from short-term to long-term, its tax falls from \(\$163.20\) to \(\$95.20\), and the optimum improves by exactly that \(\$68.00\) with the same lots. Declare a purchase of 100 shares twelve days before the sale and the wash-sale rule disallows \(\$1{,}800\) of lot 2's loss. The optimum's tax rises to \(-\$1{,}703.40\), and the disallowed loss is deferred into the basis of the purchased shares, which becomes \(\$118\).The least-tax-first rule follows N. Moehle, M. J. Kochenderfer, S. Boyd and A. Ang, "Tax-aware portfolio construction via convex optimization", Journal of Optimization Theory and Applications 189 (2021), §3.1–3.2. The rates are the top federal rates plus the \(3.8\%\) net investment income tax, 26 U.S.C. §1(h), §1(j) and §1411; the holding-period rule is 26 U.S.C. §1222 ("held for more than 1 year"); the wash-sale rule is 26 U.S.C. §1091 with the basis adjustment in §1091(d). The rates and both rules are stated from memory and remain unverified, and Section 9 gives the full statement with its caveats. The figure of Section 9 computes all 256 subsets live and matches the 100 purchased shares against the loss lots in the order of the table, lot 2 first; the regulation's own matching rule for a sale of several lots, 26 CFR §1.1091-1(b) and (c), orders them by date of acquisition and here also selects lot 2, the earliest acquired of the three loss lots sold. The example is for the integrality of lots as the simplest knapsack (Section 3.1) and for the convex-envelope relaxation and its gap (Section 2.4). It is also for the two clauses, the holding period and the wash sale, that make the tax problem of Section 9 a MINLP rather than a QP.
R6: the LP sells lots in order until the $62,000 is raised
LP x_5 = 0.3333 tax -2,246.04
optimum lots 2, 3, 4, 5, 6 and 8 whole tax -2,131.80
gap 114.24: all of lot 5, not a third of it
The two programs below reproduce the central numbers of the six examples with nothing but numpy. The first covers R1 to R3. It introduces the one tool both programs use, a vertex enumerator for a polygon given by inequalities, which is a correct linear-programming solver in two dimensions because a linear function is largest at a vertex.
# Running examples R1, R2 and R3: the central numbers the text quotes.
#
# From first principles, with numpy only. The one tool is a vertex
# enumerator for a polygon given by inequalities, a correct LP solver in
# the plane because a linear function is largest at a vertex.
import itertools
import numpy as np
def vertices(rows):
"""Vertices of {p : a.p <= b} in the plane.
The pairwise intersections of the rows that satisfy every row.
"""
V = []
for (a1, b1), (a2, b2) in itertools.combinations(rows, 2):
A = np.array([a1, a2], float)
if abs(np.linalg.det(A)) < 1e-12:
continue
p = np.linalg.solve(A, [b1, b2])
if all(np.dot(a, p) <= b + 1e-9 for a, b in rows):
V.append(p)
return np.array(V)
def lp2(c, rows):
"""Maximize c.p over the polygon: a linear function is largest at
a vertex."""
V = vertices(rows)
k = int(np.argmax(V @ c))
return V[k] @ c, V[k]
# R1: the two-variable integer program of the plane figures,
# objective c = (cos 45 deg, sin 45 deg)
P = [((2, 5), 24.5), ((5, 2), 30.5), ((-3, 4), 11), ((1, -2), 4.2),
((-1, 0), 0), ((0, -1), 0)]
c = np.array([np.cos(np.pi/4), np.sin(np.pi/4)])
zlp, plp = lp2(c, P)
pts = np.array([(x, y) for x in range(8) for y in range(7)
if all(np.dot(a, (x, y)) <= b + 1e-9 for a, b in P)])
vals = pts @ c
zint = vals.max()
# the integer points tied for best
best = [tuple(q) for q, v in zip(pts, vals) if v > zint - 1e-9]
print(f"R1: {len(pts)} integer points; "
f"LP {zlp:.3f} at ({plp[0]:.3f}, {plp[1]:.3f}); "
f"integer {zint:.3f} at {best}; gap {(zlp - zint)/zlp:.1%}")
def hull_edges(Q):
"""Andrew's monotone chain: the number of edges of the convex hull."""
Q = sorted(map(tuple, Q))
def turn(o, a, b):
return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])
def chain(S):
H = []
for p in S:
while len(H) >= 2 and turn(H[-2], H[-1], p) <= 0:
H.pop()
H.append(p)
return H[:-1]
return len(chain(Q) + chain(Q[::-1]))
print(f"R1: the convex hull of the integer points has "
f"{hull_edges(pts)} facets")
# R2: maximize x*y on the unit square under 2x + y <= 1.2; the maximum
# lies on the line y = 1.2 - 2x
xs = 1.2/4 # d/dx [x(1.2 - 2x)] = 1.2 - 4x = 0
print(f"R2: true maximum {xs*(1.2 - 2*xs):.2f} "
f"at ({xs:.1f}, {1.2 - 2*xs:.1f})")
S = [((2, 1), 1.2), ((1, 0), 1), ((0, 1), 1), ((-1, 0), 0), ((0, -1), 0)]
# McCormick on the unit square: 0 <= w, x + y - 1 <= w, w <= x, w <= y;
# so max w = max min(x, y), one LP on each side of the crease x = y
z1, p1 = lp2(np.array([1., 0.]), S + [((1, -1), 0)]) # x <= y: max x
z2, p2 = lp2(np.array([0., 1.]), S + [((-1, 1), 0)]) # y <= x: max y
zmc, pmc = max([(z1, p1), (z2, p2)], key=lambda t: t[0])
print(f"R2: McCormick relaxation {zmc:.3f} "
f"at ({pmc[0]:.1f}, {pmc[1]:.1f}), "
f"where x*y is only {pmc[0]*pmc[1]:.2f}")
def mc_max(k):
"""The same relaxation on each of the k x k sub-boxes: the largest
relaxed maximum."""
best = -np.inf
for i, j in itertools.product(range(k), repeat=2):
# the sub-box [a, b] x [lo, hi]
a, b, lo, hi = i/k, (i+1)/k, j/k, (j+1)/k
box = S + [((1, 0), b), ((-1, 0), -a),
((0, 1), hi), ((0, -1), -lo)]
# over-estimator min(lo*x + b*y - b*lo, hi*x + a*y - a*hi):
# one LP on each side of the crease
pieces = (((lo, b), -b*lo, ((lo - hi, b - a), b*lo - a*hi)),
((hi, a), -a*hi, ((hi - lo, a - b), a*hi - b*lo)))
for obj, const, crease in pieces:
V = vertices(box + [crease])
if len(V):
best = max(best, np.max(V @ np.array(obj, float)) + const)
return best
print("R2: McCormick maximum with k pieces per axis, k = 1 to 6: "
+ ", ".join(f"{mc_max(k):.3f}" for k in range(1, 7)))
# R3: st_e13, minimize 2x + y s.t. 1.25 - x^2 - y <= 0, x + y <= 1.6,
# 0 <= x <= 1.6, y in {0, 1}
for y in (0, 1):
# the smallest x the parabola allows on the row; x + y <= 1.6 is slack
x = np.sqrt(1.25 - y)
print(f"R3: y = {y}: x >= {x:.6f}, value {2*x + y:.6f}")
E = [((-1.6, -1), -1.25), ((1, 1), 1.6), ((1, 0), 1.6), ((-1, 0), 0),
((0, 1), 1), ((0, -1), 0)]
# the chord row is y >= 1.25 - 1.6x
zr, pr = lp2(np.array([-2., -1.]), E)
print(f"R3: chord relaxation at the root: "
f"bound {-zr:.4f} at ({pr[0]:.5f}, {pr[1]:.0f})")
For R1 it prints \(22\) integer points, the LP value \(5.556\) at \((4.929, 2.929)\), the integer optimum \(4.950\) attained at both \((4, 3)\) and \((5, 2)\), a gap of \(10.9\%\) and a hull with \(7\) facets. For R2 it prints the true maximum \(0.18\) at \((0.3, 0.6)\) and the McCormick value \(0.400\) at \((0.4, 0.4)\), where the product is \(0.16\). It then prints the ladder \(0.400, 0.233, 0.193, 0.192, 0.187, 0.185\) of relaxed maxima for \(k = 1, \dots, 6\) pieces per axis, each the largest value over \(k^2\) small linear programs. For R3 it prints the two rows' values, \(2.236068\) at \(x = 1.118034\) and \(2.000000\) at \(x = 0.5\), and the root bound \(1.3125\) at \((0.15625, 1)\). The cost of the vertex enumerator is cubic in the number of rows: one candidate per pair of rows, each checked against every row. The three examples are independent of one another, and the candidate vertices are independent across pairs of rows, so both would run as one batch. That is nothing here, and it is the reason the figures can afford to re-solve on every slider movement. A real LP solver is Section 7.2's subject. The second program covers R4 to R6.
# Running examples R4, R5 and R6: the central numbers the text quotes.
#
# With numpy only. vertices and lp2 are repeated from the block above so
# that this block runs on its own.
import itertools
import numpy as np
def vertices(rows):
V = []
for (a1, b1), (a2, b2) in itertools.combinations(rows, 2):
A = np.array([a1, a2], float)
if abs(np.linalg.det(A)) < 1e-12:
continue
p = np.linalg.solve(A, [b1, b2])
if all(np.dot(a, p) <= b + 1e-9 for a, b in rows):
V.append(p)
return np.array(V)
def lp2(c, rows):
V = vertices(rows)
k = int(np.argmax(V @ c))
return V[k] @ c, V[k]
# R4: Haverly's pooling problem, instance 1. With the pool quality q
# fixed the problem is linear in the flows and splits into one
# two-variable LP per product (pool flow y, direct flow z from C).
def profit(q, cB=16.0, dX=100.0):
# pool material at sulphur q: a mix of A (3%) and B (1%)
cpool = 6.0*(q - 1)/2 + cB*(3 - q)/2
zX, _ = lp2(np.array([9 - cpool, 9 - 10.0]),
[((q - 2.5, 2 - 2.5), 0), ((1, 1), dX),
((-1, 0), 0), ((0, -1), 0)])
zY, _ = lp2(np.array([15 - cpool, 15 - 10.0]),
[((q - 1.5, 2 - 1.5), 0), ((1, 1), 200.0),
((-1, 0), 0), ((0, -1), 0)])
return zX + zY
qs = np.round(np.arange(1.0, 3.0001, 0.01), 2)
prof = np.array([profit(q) for q in qs])
i = int(np.argmax(prof))
j = int(np.argmax(np.where(qs >= 2.4, prof, -np.inf)))
print("R4: profit at q = 1, 1.5, 2, 2.5, 3: "
+ ", ".join(f"{profit(q):.0f}" for q in (1, 1.5, 2, 2.5, 3)))
print(f"R4: global maximum {prof[i]:.0f} at q = {qs[i]:.2f}; "
f"the other hill tops at {prof[j]:.0f} at q = {qs[j]:.2f}")
# the McCormick relaxation's point (fA, fB, yX, yY, zX, zY, p, wX, wY):
# check its balances, its eight planes and its profit
fA, fB, yX, yY, zX, zY, p, wX, wY = 50, 100, 50, 100, 50, 100, 2.0, 150, 100
ok = [fA + fB == yX + yY, 3*fA + fB == wX + wY,
wX + 2*zX <= 2.5*(yX + zX), wY + 2*zY <= 1.5*(yY + zY),
yX + zX <= 100, yY + zY <= 200]
# McCormick planes for w = p*y on [1, 3] x [0, Y]
for w, y, Y in ((wX, yX, 100), (wY, yY, 200)):
ok += [w >= 1*y, w >= 3*y + Y*p - 3*Y, w <= 3*y, w <= 1*y + Y*p - 1*Y]
print(f"R4: McCormick point feasible for all {len(ok)} rows: {all(ok)}; "
f"profit {9*(yX+zX) + 15*(yY+zY) - 6*fA - 16*fB - 10*(zX+zY):.0f}; "
f"it carries sulphur {wX/yX:.0f}% to X and {wY/yY:.0f}% to Y "
f"from one pool")
# R5: three assets, at most two held: the least variance that reaches a
# target return R, with and without the limit
MU = np.array([0.10, 0.07, 0.04])
SIG = np.array([0.20, 0.15, 0.10])
rho, R = 0.3, 0.07
S = (1 - rho)*np.diag(SIG**2) + rho*np.outer(SIG, SIG)
def minvar(support):
"""min w'Sw s.t. 1'w = 1, mu'w >= R, w >= 0 on a support, by every
active set."""
best, idx = np.inf, list(support)
k = len(idx)
Ss = S[np.ix_(idx, idx)]
for ret in (False, True):
for r in range(k):
for zeros in itertools.combinations(range(k), r):
A = np.array([np.ones(k)] + ([MU[idx]] if ret else [])
+ [np.eye(k)[z] for z in zeros])
m = len(A)
rhs = [1.0] + ([R] if ret else []) + [0.0]*len(zeros)
K = np.block([[2*Ss, A.T], [A, np.zeros((m, m))]])
try:
w = np.linalg.solve(K, np.r_[np.zeros(k), rhs])[:k]
except np.linalg.LinAlgError:
continue
if (w >= -1e-9).all() and MU[idx] @ w >= R - 1e-9:
best = min(best, w @ Ss @ w)
return best
v_all = minvar((0, 1, 2))
v_two = min(minvar(s) for s in itertools.combinations(range(3), 2))
print(f"R5: at R = {R:.0%}: no limit on the number of assets, "
f"sigma = {np.sqrt(v_all):.2%}; "
f"at most two assets, sigma = {np.sqrt(v_two):.2%}")
# R6: eight tax lots of one stock at $100 a share, $62,000 to raise:
# every subset of whole lots, then the LP relaxation
sh = np.array([100, 150, 80, 120, 60, 200, 90, 50])
basis = np.array([62, 118, 95, 131, 88, 104, 71, 99])
days = np.array([900, 400, 340, 45, 500, 200, 1200, 30])
price, cash = 100.0, 62000.0
# long-term after one year: 23.8% against 40.8%
rate = np.where(days > 365, 0.238, 0.408)
tax, proceeds = sh*(price - basis)*rate, sh*price
best, feasible = None, 0
for m in range(256):
L = [i for i in range(8) if m >> i & 1]
if proceeds[L].sum() >= cash:
feasible += 1
if best is None or tax[L].sum() < best[0]:
best = (tax[L].sum(), L, proceeds[L].sum())
print(f"R6: {feasible} of 256 subsets raise enough; "
f"optimum tax {best[0]:,.2f} "
f"selling lots {[i+1 for i in best[1]]}, raising ${best[2]:,.0f}")
# the LP: loss lots whole, then gain lots in increasing order of tax per
# dollar, the last in part
x = (tax <= 0).astype(float)
pr, tx = (x*proceeds).sum(), (x*tax).sum()
for i in sorted(np.where(tax > 0)[0], key=lambda i: tax[i]/proceeds[i]):
if pr >= cash:
break
f = min(1.0, (cash - pr)/proceeds[i])
x[i] = f
pr += f*proceeds[i]
tx += f*tax[i]
print(f"R6: LP relaxation tax {tx:,.2f} raising ${pr:,.0f}, "
f"lot 5 sold in part (x_5 = {x[4]:.4f}); "
f"gap {best[0] - tx:,.2f}")
For R4 it prints the profits \(400, 300, 0, 50, 100\) at \(q = 1, 1.5, 2, 2.5, 3\), the global maximum \(400\) at \(q = 1.00\) and the other hill's top \(100\) at \(q = 3.00\). It also checks that the McCormick relaxation's point satisfies all fourteen rows with profit \(500\) while carrying \(3\%\) sulphur to X and \(1\%\) to Y from one pool. That \(500\) is the relaxation's optimum, and not merely a value it allows, is a nine-variable linear program, which the figure of Section 3.5 solves. For R5 it prints the volatilities \(11.12\%\) without the limit and \(12.45\%\) with it. The perspective bound \(11.71\%\) is a small conic program and is quoted from the figure of Section 4.4. For R6 it prints \(34\) feasible subsets, the optimum \(-2{,}131.80\) selling lots 2, 3, 4, 5, 6 and 8 for \(\$66{,}000\), and the relaxation \(-2{,}246.04\) with \(x_5 = 0.3333\) and a gap of \(114.24\). Each example is solved here by brute force: a sweep over \(201\) values of \(q\), an enumeration of active sets, an enumeration of \(256\) subsets. The sweep and the enumerations are independent evaluations that a GPU would run as one batch. Sections 2 to 4 replace such enumerations, which grow as \(k^d\) or \(2^p\), by bounds that discard most candidates without evaluating them.
Where this is used
The six examples stand for six structures a solver meets. R1 is a lattice inside a polytope. R2 is a product of two variables. R3 is a constraint that bends the wrong way with a binary variable beside it. R4 is a network of bilinear balances. R5 is a set of indicators attached to a convex quadratic, and R6 is a knapsack with rules attached. Every global solver named in Section 5 handles the first by the machinery of Sections 3.1 to 3.3, the second to fourth by the envelopes and spatial branching of Sections 2.4 and 3.5, and the fifth by the perspective and conic reformulations of Section 4. The sixth is the shape of the problem Section 9 sets out to solve.
What parallelizes
For each of the six, the computation that dominates, and the one a GPU is asked to accelerate, is the repeated solution of a convex relaxation over many boxes or many subproblems at once. For R2 these are the sub-boxes, for R3 the one-dimensional nodes, for R4 the quality intervals, for R5 the supports and for R6 the sub-knapsacks. Sections 6 and 7 are about doing that in batches, and about the price, in the validity of the bounds, of doing it with the inexact arithmetic a device offers.